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ABSTRACT 


A  framework  for  modelling  and  simulation  of  phenomena 
in  developmental  genetics  is  presented  in  this  thesis. 
Initially,  formalisms  for  such  modelling  are  introduced, 
using  systems  theory  to  provide  a  context  for  questions 
concerning  model  validity. 

The  biological  and  theoretical  background  of  the 
project  is  discussed.  Phenomena  such  as  transdetermination, 
compar tmenta 1 ization  and  positional  information  are 
described.  Mechanisms  for  pattern  formation  such  as 
chemical  reaction  and  diffusion  systems  are  investigated. 
Existing  computer -or i ented  models  are  surveyed.  The 

experimental  frame  for  the  project  is  discussed. 

An  underlying  model  for  development  is  presented,  and  a 
modelling  framework  established,  so  that  particular 
phenomena  can  be  modelled  by  selecting  options  from  a  model 
pool.  Structures  are  defined  to  enable  representat i on  of 
individual  cells,  groups  of  cells  in  structural 
relationships,  and  developmental  systems.  The  modelling 
framework  is  mu  1 1 i 1  eve  1 1 ed  and  of  mixed  type. 

The  implementation  of  the  modelling  framework  in  ALGOLW 
is  described.  The  validity  of  approximations  necessary  in 
representing  the  modelling  framework  in  iterative  form  is 
discussed.  Such  approximations  include  both  discretization, 
of  the  time  base  and  continuous  state  variables,  and 
replacement  of  parallel  by  sequential  processing. 
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Particular  experiments  performed  in  developing 
modelling  options  are  described.  The  control  of  growth  is 
investigated,  and  options  provided  for  locally  mediated  cell 
division.  Intercel lular  communication  by  means  of  reaction 
and  diffusion  is  modelled  by  means  of  numerical 
approximation,  and  it  is  shown  that  the  concentration 
patterns  predicted  theoretically  are  observed  at  critical 
sizes  of  cell  structures.  The  option  of  modelling  the 
sequence  of  patterns  by  means  of  analytic  solutions  is 
provided  for  rectangular  structures.  A  model  is  given  for 
determination,  wherein  a  cell's  genetic  information  is 
affected  on  recognition  of  a  critical  value  in  its  local 
conditions.  Finally,  models  for  segmentation  in  the 
blastoderm  of  Drosophila  and  for  the  development  of  the 
chick  wing  are  presented,  using  options  within  the  modelling 
framework . 

Contributions  to  developmental  genetics  and  to  the  art 
of  modelling  and  simulation  are  discussed.  Finally,  the 
formal  context  employed  in  developing  the  modelling 
framework  is  evaluated  in  terms  of  its  success  in  providing 
a  mathematical  environment  for  the  study  of  developmental 
genet i cs . 
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CHAPTER  ONE 

Introduct ion 


1 . 1  Mot i vat i on 

Developmental  genetics  is  a  vital  area  in  current 
biological  research.  The  genetic  code  has  been  deciphered, 
and  investigations  proceed  into  the  means  by  which  this 
inherited  information  stored  in  the  cell  is  translated  into 
biochemical  activity.  In  a  developing  organism,  this 
translation  of  genetic  information  results  in  a  sequence  of 
changes  in  cells  and  their  progeny  leading  to  a  cooperative 
population  of  many  cell  types,  that  is,  a  viable  organism. 
The  genetic  material  in  a  cell  may  be  regarded  as  a 
specification  of  all  chemical  activities  possible  in  this 
organism:  during  development,  each  cell  has  operative  only 
parts  of  this  material.  Which  part  of  the  genetic  material 
is  operative  at  a  given  time  depends  on  the  operational 
history  of  the  cell  (its  initial  state  and  all  succeeding 
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states),  any  external  influences,  and  its  position  and 
relationship  to  other  cells  in  the  organism.  An  immediate 
comparison  may  be  made  to  the  operation  of  an  algorithm  in  a 
computer:  which  portion  of  the  algorithm  is  in  action 
depends  on  operation  of  the  algorithm  to  date  and  inputs  to 
which  the  computer  may  be  connected. 

The  analogy  between  developing  organism  and  computer  is 
developed  into  a  simulation  project  in  this  thesis.  The 
project  is  a  response  to  certain  needs  both  in  genetics  and 
in  that  area  of  computing  science  Known  as  modelling  and 
simulation. 

1.1.1  The  Role  of  Mathematics  in  Genetics 

In  the  physical  sciences  mathematics  has  played  an 
essential  role  in  providing  a  formal  and  exact  context  for 
the  formulation  of  theories  and  their  experimental 
verification.  The  aim  of  biological  science  is  similar  to 
that  of  the  physical  sciences,  that  is,  the  development  of 
generalized,  unifying  theories  about  its  own  domain.  While 
much  of  modern  mathematics  has  developed  in  interaction  with 
the  physical  sciences,  many  important  milestones  in  biology, 
and  in  particular  in  genetics,  have  been  achieved  when 
conventional  mathematics  has  been  applied  in  biology. 

Scientific  biology  may  be  considered  to  have  begun  when 
it  was  recognized  that  the  metabolic  functions  and 
transformat i ons  in  organisms  are  reducible  to  simple 
chemistry.  The  modern  biochemistry  of  enzymatic  reactions, 
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based  on 

the  work  of 

Ha  1 dane 

(  1930 

) ,  cons i sts 

of 

quant i tat i ve 

reaction  laws 

based  on 

the 

mathemat i cs 

of 

physical  chemistry.  Differential  equations  express  changes 
of  concentration  of  interdependent  chemical  species.  Thus, 
at  the  molecular  level,  biology  has  a  useful  mathematical 
forma  1 i sm . 

The  mathematical  theory  of  evolution  and  natural 
selection,  in  which  the  pioneers  were  Wright  (1931),  Haldane 
(1932)  and  Fisher  (1930),  rests  on  the  recognition  of  the 
gene  as  a  heritable  unit  and  the  application  of  stochastic 
mathematics  to  genetic  populations.  Thus  at  the  level  of 
the  genetic  population  there  is  a  useful  mathematical 
forma  1 i sm . 

Consideration  of  possible  mechanisms  by  which  genetic 
information  may  be  interpreted  differentially  led  to  the 
discovery  of  phenomena  such  as  end  product  inhibition,  in 
which  the  end  product  of  a  synthetic  reaction  combines  with 
a  catalyst  of  the  same  reaction  to  modify  the  reaction  rate. 
Thus  control  systems  theory  began  to  be  used  in  the  context 
of  genetic  mechanisms,  and  in  fact  attempts  have  been  made 
to  describe  the  whole  organism  as  an  adaptive  control  system 
(Reiner,  1968;  Rosen,  1972).  Thus  both  at  the  molecular  and 
at  the  organic  level ,  systems  theory  provides  a  useful 
mathematical  formalism.  The  use  of  this  formalism  will  be 
discussed  in  greater  detail  in  Chapter  Two. 

Dynamic  processes  by  which  form  and  pattern  are 
generated  in  biological  systems  may  be  grouped  under  the 
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name  morphogenesis.  The  fundamental  question  in  such 
processes  is  how  patterns  arise  autonomously  in  initially 
uniform  situations.  A  theory  first  developed  by  Rashevsky 
(1960)  and  Turing  (1952)  suggests  that  the  uniform 
configuration  must  be  the  result  of  an  unstable  equilibrium 
which  may  be  transformed  by  random  fluctuations  into  a 
nonuniform  stable  state.  Mathematical  formalisms  which 
attempt  to  deal  with  such  phenomena  include  bifurcation 
theory  (Nicolis  and  Prigogine,  1977),  catastrophe  theory 
(Thom,  1972),  and  stability  theory  (Gmitro  and  Scriven, 
1966).  These  theories  will  be  discussed  in  greater  detail 
in  Chapter  Three.  In  developmental  genetics  the  scope  of 
such  a  treatment  will  usually  be  the  morphogenetic  field, 
which  is  of  the  order  of  100  cells  (Wolpert,  1978;  Crick, 
1970) . 

There  are  many  examples  of  attempts  to  apply 
quantitative  methods  to  individual  phenomena  in 
developmental  genetics;  some  of  these  will  be  discussed  in 
Chapter  Three. 

It  appears  likely,  therefore,  that  any  single 
phenomenon  in  genetics  may  be  explicable  in  terms  familiar 
to  the  mathematics  of  physical  science.  At  the  molecular 
level  there  is  enzyme  kinetics;  at  the  level  of  the  cell 
population  there  is  statistics;  at  the  level  of  the 
morphogenetic  field  there  may  be  a  mathematical  formalism 
for  equilibrium  states.  However  development  is  a  phenomenon 
occurring  at  all  these  levels,  and  it  has  an  additional 
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characteristic.  It  is  a  process  occurring  in  and  affecting 
a  structure  of  cells. 

Essential  features  of  pattern  formation  in  a  growing 
organism  are  the  relative  positions  of  cells  and  their 
intercommunication.  The  mathematics  of  classical  mechanics 
deals  with  points,  while  that  of  modern  physics  deals  with 
fields;  neither  context  is  sufficient  to  describe 
development . 

The  computing  scientist  may  suggest  that  the  "cellular 
space",  wherein  independently  functioning  "cells"  are  tied 
into  a  structure  by  means  of  connections  to  neighbours,  is 
the  mathematical  tool  suited  to  development.  In  addition 
there  is  an  extant  body  of  mathematics  related  to  cellular 
automata;  however  there  are  features  of  a  developing 
organism  which  make  conventional  cellular  automata  theory 
insufficient.  Prominent  among  these  are,  firstly, 

nonuniformity  of  function  in  the  cells  of  increasing 
complexity  as  cells  differentiate,  local  functioning  not 
being  subject  to  local  control,  and  secondly,  the  continual 
growth  and  change  in  the  structure  as  cells  divide  or  move. 
The  structure  of  cells  is  dynamically  subject  to  controls 
which  are  mediated  at  the  level  of  the  organism,  or  at  least 
the  morphogenetic  field. 

A  mathematical  context  for  development  must  be  capable 
of  mirroring  the  interaction  between  process  (the  dynamic 
aspect)  and  structure  (the  static  aspect).  Furthermore,  it 
should  be  capable  of  invoking  the  mathematical  tools 
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mirroring  events  at  molecular  and  collective  levels.  This 
thesis  outlines  the  development  of  a  project  in  which 
systems  theory  is  used  in  alliance  with  a  cellular  space  to 
attempt  this  task. 

An  individual  biological  cell  can  act  as  an  automaton 
(for  example,  in  determining  and  retaining  a  cell  type 
marker )  or  as  a  location  in  a  field  (for  example,  in 
determining  local  concentration  of  a  freely  flowing 
chemical).  Thus  a  developing  organism  has  both  continuous 
and  discrete  aspects  (Lewis,  1976).  The  operation  of  field¬ 
like  phenomena  can  "feed  back"  to  change  structure  and  local 
functioning  in  the  cellular  space.  While  the  technique  will 
not  be  discussed  in  greater  detail  here,  two  observations 
may  be  made.  Firstly,  a  heirarchy  of  control  is  implied, 
that  is,  the  interaction  of  process  and  structure  occurs  at 
a  non-local  or  "higher"  level.  Secondly,  it  is  likely  that 
various  options  should  be  present  in  any  situation,  for 
example,  different  mathematical  theories  of  cell  division  or 
of  cell  intercommunication  may  be  desirable.  A  "pool"  of 
techniques  should  therefore  be  provided. 

A  computer  is  an  ideal  tool  in  a  situation  where  many 
cells,  each  capable  of  various  behaviors,  interact  in  a 
dynamic  fashion.  Its  large  storage  capacity  allows 
consideration  of  many  factors,  while  its  capacity  to 
generate  data  dynamically  following  an  iterative  program 
specification  allows  for  dynamic  testing  of  theories  of 
development . 
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Developmental  genetics,  in  summary 
the  evolution  of  mathematical  formalisms 
particular  features.  When  individual 
related  to  such  a  context,  more  mean 
theories  of  development  will  result, 
field  where  the  number  of  parameters  and 
more  fruitful  direction  for  experimentat 
in  a  unifying  context  (Reiner,  1968 
establishes  a  mathematical  context 
experiments  which  includes  computer  simul 
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options  is  vast, 
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for  developmental 
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1.1.2  Scientific  Modelling  and  Simulation 

A  hypothesis  about  a  natural  phenomenon  may  be 
expressed  in  mathematical  (formal,  symbolic,  and  exact) 
terms.  How  is  this  formulation  related  to,  firstly,  the 
phenomenon  it  is  supposed  to  represent,  and  secondly,  the 


refinements  through  which  the 

testing  the  hypothesis? 

exper i menter 

may 

put  it  in 

A  natural  phenomenon 

is  Known  to 

the 

exper imenter 

solely  through  its  behavior 

( measurements 

on 

parameters , 

response  to  stimuli  or 

inputs,  etc 

)  . 

He  makes  a 

hypothetical  construct  which,  he  hopes,  represents  reality 
in  some  fashion;  usually,  it  is  designed  to  generate  similar 
behavior . 

Modelling  and  simulation  are  terms  used  in  describing 
the  process  wherein  a  hypothesis  is  translated  into  formal 
terms  which  may  then  be  used  in  an  iterative  manner  to 
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generate  artificial  behavior  (to  be  compared  to  the  real 
data  available  to  the  exper i menter ) .  The  process  becomes  a 
science  when  the  relationship  between  model  and  reality,  and 
between  one  modelling  formalism  and  another,  can  be  given  in 
equally  formal  terms.  Only  thus  can  the  intended  scope  of  a 
model  be  defined  and  its  significance  properly  assessed. 

Ziegler  (1976)  has  proposed  a  theory  of  modelling  and 
simulation  which  suggests  the  manner  in  which  a  modelling 
experiment  should  be  constructed  and  formalized  if  questions 
of  validity  are  to  be  addressed  directly.  Modelling  in 
developmental  genetics,  with  its  multiplicity  of  parameters 
and  influences,  appears  an  especially  challenging  task. 
This  thesis  will  use  Zeigler's  methods  wherever  possible, 
both  as  an  aid  in  a  complex  project  and  in  an  attempt  to 
assess  the  applicability  of  the  theories.  Zeigler's 
theories  are  particularly  suited  to  this  project  because 
they  are  couched  in  systems  theoretic  terms;  thus  the 
dynamic  aspects  of  the  project  can  be  included  naturally. 

1.1.3  Relations  to  Computing  Science 

Simulation  is  a  technique  used  widely  by  computing 
scientists,  whether  in  the  context  of  efficient  data 
processing  (for  example,  queueing  models)  or  in  more  applied 
fields  (for  example,  information  retrieval  or  ecological 
models).  A  project  such  as  this  one,  though  specialized  in 
aim,  will  provide  contributions  to  the  general  art  of 
simulation  itself. 
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The  project  is  concerned  with  dynamically  occurring 
changes  in  structure  and  local  function  in  a  cell  space.  In 
simulation  and  in  other  areas  of  computing  science  where  the 
cell  space  formalism  is  used,  considerable  power  and 
extension  of  applicability  can  be  expected  if  methods  of 
dealing  with  change  and  non-uniformity  are  developed. 

Biological  information  processing  systems  provide  a 
rich  environment  for  both  mu  1 1 i 1  eve  1 1 ed  and  mixed- type 
modelling.  In  a  developmental  context  modelling  can  be  done 
at  the  various  levels  of  the  organism,  part  of  the  organism, 
the  morphogenetic  field,  the  cell,  or  the  i ntr ace  1 1 u 1 ar 
genetic  apparatus,  or  at  a  combination  of  these  levels.  At 
different  levels  different  modelling  techniques  may  be 
appropriate;  furthermore,  it  may  be  possible  to  model  the 
same  phenomenon  at  varying  levels  of  "resolution",  with 
different  techniques,  and  thus  have  a  means  of  checking 
results  (Sampson,  in  prepar at i on ) .  As  discussed  above,  a 
developmental  system  contains  both  continuous  elements  (for 
example,  time  and  field-like  influences)  and  discrete 
elements  (for  example,  certain  cellular  char acter i s t i cs ) , 
and  thus  must  be  modelled  as  a  system  of  mixed  type. 
Paradigms  for  mu  1 1 i 1  eve  1 1 ed ,  mixed  type  modelling  are 
essential  for  many  real  systems  of  interest,  and  therefore 
successes  in  these  areas  may  be  expected  to  contribute  to 
the  art  of  simulation. 

Computing  science  is  the  science  of  successful 
information  processing.  Biological  information  processing 
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has  been  proven  successful  through  the  process  of  evolution, 
and  therefore  may  have  ideas  to  contribute  to  artificial 
systems.  In  particular,  biological  systems  store 
information  compactly  (the  memory  element  being  at  the 
molecular  level),  are  reliable  for  the  lifespan  of  the 
system,  and  above  all  behave  adaptively.  Developmental 
systems  exhibit  particularly  powerful  adaptation  on  the 
basis  of  experience,  as  cells  and  structures  develop 
specialized  functions.  Artificial  information  processing 
looks  toward  more  compactness,  greater  reliability,  and 
adaptibility  in  both  hardware  and  software.  While  the 
mechanisms  of  biological  systems  may  not  be  directly  and 
currently  applicable  in  artificial  systems,  an  understanding 
of  such  mechanisms  may  be  expected  to  contribute  to  the 
design  of  the  artificial  system  of  the  future. 

In  summary,  then,  simulation  in  developmental  genetics 
has  special  features  and  difficulties  which,  if  overcome, 
may  be  expected  to  contribute  both  to  the  art  of  simulation 
and  to  computing  science. 

1.2  Organization  of  the  Thesis 

Since  Zeigler' s  theory  of  modelling  and  simulation 
(Zeigler,  1976)  is  used  as  a  formal  background  in  developing 
this  project,  Chapter  Two  is  devoted  to  a  concise  exposition 
of  the  theory.  The  theory  divides  such  an  enterprise  into 
five  elements:  the  real  system,  the  experimental  frame,  the 
base  model,  the  lumped  model,  and  the  computer.  Chapters 
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Three,  Four  and  Five  explore  these  elements  in  the  context 
of  the  present  project. 

Chapter  Three,  entitled  The  Real  System,  introduces  the 
necessary  genetics  and  the  actual  areas  in  which  the  project 
operates.  Since  the  state  of  Knowledge  and  direction  of 
experimentation  in  the  subject  seems  inextricable  from 
models  under  current  consideration,  the  nature  of  existing 
models  is  also  described  in  this  chapter;  these  include  both 
purely  genetic  modelling  and  work  involving  simulation.  The 
chapter  ends  by  discussing  means  by  which  experimental 
frames  will  be  selected  for  particular  modelling 
exper i ments . 

Chapter  Four  describes  the  base  model  for  developmental 
systems  (the  complete  hypothetical  model  capable  of 
generating  all  possible  behaviors  in  all  experimental 
frames).  The  modelling  framework  is  described  in  informal 
and  formal  terms.  Those  elements  of  particular  lumped 
models  (simplified  versions  of  the  base  model  viewed  in 
particular  experimental  frames)  which  are  held  in  common  are 
then  described  informally  and  formally. 

In  Chapter  Five  translation  of  the  lumped  models  into 
iterative  form  is  discussed.  Such  issues  as  the  choice  of 
programming  environment,  treatment  of  the  time  base,  data 
structures,  and  problems  of  economization  are  discussed. 

Chapter  Six  describes  the  implementation  of  particular 
lumped  models  (models  of  certain  developing  systems).  The 
results  of  simulations  are  presented  and  discussed,  with 
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attention  to  questions  of  validity  and  general  implications. 

Chapter  Seven  evaluates  the  success  with  which  the 
various  aims  of  the  project  have  been  achieved. 

1.3  Summary  of  Goals 

The  most  immediate  aim  of  the  project  is  to  provide  a 
useful,  functioning  framework  for  testing  models  in 
developmental  genetics.  Achievement  of  this  goal  will 
depend  on  reaching  many  small  individual  goals  in  obtaining 
successful  results  in  particular  experiments.  It  is  only  in 
the  context  of  such  success  that  larger  implications  of  the 
project  may  be  discussed.  These  will  include  evaluation  of 
the  extent  to  which  a  satisfying  mathematical  context  for 
development  has  been  achieved,  and  demonstration  of  the  role 
which  the  computer  can  play  in  such  a  context.  In  addition 
it  wi 1 1  be  possible  to  assess  Zeigler's  theory  of  modelling 
and  simulation  as  a  practical  strategy  in  projects  of  this 
complexity.  Further,  contributions  to  the  art  of  simulation 
and  to  computing  in  general  are  expected. 
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CHAPTER  TWO 

Modelling  and  Simulation 

Three  components  of  the  process  Known  as  "modelling  and 
simulation"  are  the  real  system,  the  model,  and  the 
computer.  A  model  is  constructed  as  a  conjecture  about  the 
nature  of  a  real  system,  and  a  simulation  built  to  generate 
the  behavior  of  the  model.  Zeigler  (1976)  has  developed  a 
framework  for  defining  these  components  and  establishing 
relationships  between  them  in  theoretical  terms.  The  essence 
of  his  treatment  is  presented  here,  with  the  aim  of 
introducing  the  present  author's  models  and  simulation 
experiments  in  perspective. 

2.1  The  Five  Elements 

Zeigler  expands  the  three  components  to  five 
conceptually  distinct  elements  of  modelling  and  simulation. 
These  are  the  real  system,  the  experimental  frame,  the  base 
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model,  the  lumped  model,  and  the  computer. 

2.1.1  The  Rea  1  System 

The  real  system  is  a  natural  or  artificial  source  of 
behavioral  data.  The  data  are  specified  as  values  of 
descriptive  variables,  which  may  be  observable  (measurable 
by  current  techniques)  or  nonobservable  (inaccessible  to  the 
i nvest i gator ) . 

Observable  variables  are  further  classified  as  either 
input  variables  (those  getting  their  values  from  influences 
external  to  the  system)  or  output  variables  (those  receiving 
values  as  the  result  of  values  of  input  variables).  Over  a 
given  period  of  time,  sequences  of  values  of  all  observable 
variables  may  be  collected  into  input  and  output  segments 
for  each  variable.  An  I/O  segment  pair  is  the  combination  of 
an  input  segment  and  an  output  segment  for  a  certain  time 
interval.  The  set  of  all  possible  I/O  segment  pairs  is  all 
that  can  be  directly  Known  about  the  real  system,  and  is 
Known  as  the  "I/O  behavior"  of  the  system. 

2.1.2  The  Experimental  Frame 

Complexity  of  the  real  system  or  limited  resources 
available  to  the  experimenter  usually  limit  access  to  the 
real  system  to  be  studied.  Associated  with  this  limiting 
experimental  frame  is  a  subset  of  the  I/O  behavior  of  the 
real  system. 
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2.1.3  The  Base  Model 

The  base  model  is  a  set  of  instructions  capable  of 
generating  all  the  I/O  behavior  of  the  real  system.  All 
variables  in  all  possible  experimental  frames  are  accounted 
for.  In  other  words,  the  base  model  is  a  hypothetical 
complete  explanation  of  the  real  system. 

2.1.4  The  Lumped  Model 

Within  an  experimental  frame,  models  may  be  constructed 
which  are  less  complex  than  the  base  model  and  yet  capable 
of  generating  the  subset  of  I/O  behavior  of  that  frame. 
Various  simplification  procedures  may  be  used  to  produce  a 
lumped  model  from  the  base  model;  these  include  eliminating 
descriptive  variables,  approximating  deterministic,  causal 
sequences  affected  by  many  factored  stochastic  events, 
coarsening  the  range  of  descriptive  variables,  and  pooling 
(or  "lumping")  components  and  their  descriptive  variables. 
The  hazards  of  simplification  and  degrees  of  validity  of  the 
resulting  lumped  models  will  be  discussed  when  a  formal 
terminology  for  model  description  has  been  introduced. 

2.1.5  The  Computer 

The  computer  is  a  device  which  can  take  a  suitable 
specification  of  the  lumped  model  and  generate  its  I/O 
pairs.  The  iterative  execution  of  the  model  instructions  is 
Known  as  simulation. 

In  terms  of  these  five  elements,  "modelling"  is  the 
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process  of  constructing  a  lumped  model  which  is  valid  with 
respect  to  the  real  system  in  the  context  of  an  experimental 
frame.  "Simulation"  is  the  process  of  generating  the 
behavior  of  a  model  by  means  of  the  execution  of  a  computer 
progr am . 

The  computer  program  itself  can  be  regarded  as  a 
specification  of  the  lumped  model.  However,  such 
specification  is  machine  and  language  dependent,  and  its 
structure  may  reflect  computational  facility  rather  than  the 
inherent  nature  of  the  model.  It  is  desirable  to  obtain  a 
formal  model  description  which  is  directly  imp lemen table  and 
also  provides  a  means  of  comparison  between  models  (for 
example,  when  simplification  procedures  are  being  employed 
to  produce  more  lumped  models).  Furthermore,  if  any  model 
can  be  expressed  in  a  standard  form,  its  salient 
distinguishing  features  can  be  immediately  recognized  and 
decisions  regarding  means  of  simulation  made  most 
efficiently.  Zeigler  describes  a  hierarchy  of  model 
specifications  of  increasing  rigour,  from  verbal 
descriptions  of  components  with  their  descriptive  variables 
and  interactions,  to  "iterative  system  specifications".  The 
latter  are  system  theoretic  specifications  and  amount  to 
description  of  a  model  in  a  certain  canonical  form.  The 
construction  of  such  a  specification  is  described  below. 
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2.2  Formal  System  Specification 

A  dynamic  system  is  a  set  of  components  with  attributes 
which  vary  with  time  and  interact  with  each  other.  The 
following  fundamental  elements  define  a  system. 

2.2.1  Time  Base 

A  time  base  is  a  set  T  which  is  isomorphic  either  to 
the  reals,  in  which  case  T  is  said  to  be  continuous,  or  to 
the  integers,  in  which  case  T  is  said  to  be  discrete.  In 
both  continuous  and  discrete  form,  T  respectively  obeys  the 
following  properties  of  the  reals  and  the  integers. 

Linear  ordering  allows  definition  of  past,  present  and 
future.  The  stipulation  that  T  be  an  Abelian  group  under 
addition  allows  consideration  of  translation  of  events  in 
time.  The  requirement  that  T  be  unbounded  above  and  below 
allows  consideration  of  arbitrarily  large  time  segments. 
Chronology  of  events  is  preserved  (under  translation)  if 
addition  is  required  to  be  order  preserving 
( 1 1  <  t2  -->  tl  +  t3  <  t2  +  1 3  for  all  t3  6  T). 

2.2.2  State  Set 

Associated  with  the  components  of  the  system  are 
descriptive  variables  which  fully  define  them.  At  any  point 
in  time,  the  state  of  the  system  is  given  by  the  values  of 
these  variables.  The  state  variable  set,  Q,  is  the  minimal 
set  which  specifies  the  system  as  it  undergoes  changes.  The 
values  of  these  var i ab 1 es , together  with  the  inputs,  at  any 
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time  t  determine  (uniquely  in  the  case  of  models  with  no 
stochastic  elements)  the  values  of  all  descriptive  variables 
at  all  future  t imes . 

2.2.3  Input  and  Output  Sets 

A  nonautonomous  system  exists  in  some  environment,  with 
which  it  interacts  via  inputs  and  outputs.  The  input  set,  X, 
is  a  collection  of  variables  which  influence,  but  are  not 
influenced  by,  the  system.  The  output  set,  Y,  is  some  subset 
of  the  descriptive  variables  of  particular  interest  in  the 
given  system. 

The  state,  input  and  output  variables  may  be  finite, 
discrete,  continuous  or  mixed,  and  the  time  base  may  be 
discrete  or  continuous.  Thus  a  system  may  be  categorized  by 
the  nature  of  its  time  base  and  variables. 

2.2.4  System  Behavior 

A  system  in  some  state  at  a  given  time  experiences  a 
set  of  input  values.  It  undergoes  a  transformat  ion  to  a  new 
state  and  produces  certain  output  values.  Such  behavior  can 
be  described  in  terms  of  two  functions,  the  transition 
function,  f,  and  the  output  function,  g: 

f:  Q  x  X  -->  Q 
g:  Q  -->  Y 

The  functions  f  and  g  may  be  deterministic,  if  each 
element  i n  Q  x  X  or  Q  maps  to  only  one  element  of  Q  or  Y ,  or 
nonde termi n i s t i c ,  if  some  domain  elements  each  map  to  more 
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than  one  range  element.  If  the  choice  of  range  element  is 
determined  stochastically,  then  the  system  is  termed 
stochastic . 

A  complete  system  description,  then,  is  given  by  the 
structure  <T , X , Q , Y , f , g> .  When  some  of  the  elements  are  not 
Known,  or  not  Known  completely,  a  less  detailed  structure 
must  be  used  to  describe  the  system. 

A  "model"  is  any  system  which  corresponds  in  some  way 
to  another  system.  The  modeller's  intention  is  to  develop  a 
fully  specified  system  as  a  model  of  one  which  is  not  fully 
understood  and  therefore  not  fully  specified.  For  purposes 
of  simulation,  the  complete  specification  of  a  model  should 
be  in  iterative  form,  that  is,  capable  of  step-by-step 
computation  from  an  initial  specified  state.  A  hierarchy  of 
system  specifications  can  be  described  in  which  an 
increasing  amount  is  Known  or  specified,  culminating  in 
iterative  specification. 

2.3  Manipulation  of  the  Time  Base 

Inherent  in  the  notion  of  iterative  specification  is 
the  idea  of  manipulations  on  the  time  base.  At  each 
computational  instant,  inputs  are  treated  as  arriving  now 
(t  =  0),  and  state  transitions  are  immediate.  Furthermore, 
only  the  behavior  relevant  to  the  current  time-step  length 
is  to  be  considered.  As  a  preliminary  to  the  description  of 
specification  hierarchies,  operations  on  and  assumptions 
about  the  time  base  will  be  presented. 
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Let  Z  be  one  of  X,  Q  and  Y.  Let  < 1 1  ,  1 2 >  be  a  time 
interval.  Here  <  >  implies  (  ),  (  ],[  ) ,  or  [  ].  Then  let  w 
be  a  mapping  from  <t1,t2>  to  Z,  that  is,  an  input  or  output 
segment  or  a  state  trajectory  from  tl  to  t2. 

w:  <t1 , 1 2 >  -->  Z 

w  is  Known  as  a  segment  in  (Z,T).  If  wl ,  w2  are 
segments  in  (Z,T)  and  wl:  <  1 0 , 1 1 >  -->  Z  and  w2  : 
< 1 1  , 1 2 >  -->  Z  they  are  said  to  be  contiguous.  A  composition 
operator  (•)  may  be  defined  for  contiguous  segments: 

w1»w2  =  wl  for  t  6  < 1 0  ,  1 1 > 
w2  for  t  6  <t1, t2> 

If  a  unique  value  is  assigned  at  tl  (according  to 
whether  >  means  )  or  ] ,  and  <  means  (  or  [ )  then  composition 
is  associative:  let  wl  and  w2  be  contiguous,  and  also  w2  and 
w3 .  Then 

(w1»w2)«w3  =  w1*(w2*w3) 

The  notion  of  composition  allows  formalisation  of  the 
idea  of  successive  input  experiments.  A  segmentation 
operator  may  also  be  defined,  allowing  experiments  to  be 
decomposed  into  smaller  input  setments.  In  addition,  a 
translation  operator  allows  a  segment  to  be  applied  at  a 
d i f f erent  time. 

Given  a  set  of  segments  W,  a  subset  of  (Z,T),  it  is 
desired  to  find  a  set  of  segments  G  such  that  compositions 
of  segments  in  G  yield  the  segments  of  W.  Such  a  set  G  is 
called  a  set  of  generators  for  W,  or  G  is  said  to  generate 
W.  The  set  {wl ,w2 , . . . , wn}  is  a  decomposition  of  w  by  G  if 
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wi  6  G,  i  =  1,2,  .  .  . ,n,  and  w  =  w1*w2«...wn. 

In  the  interests  of  iterative  representat i ons  it  is 
desirable  to  obtain  unique  (canonical)  decompositions  of 
segments.  Such  decompositions  may  be  obtained  by  a  process 
called  maximal  length  segmentation  (mis). 

The  process  of  finding  an  mis  decomposition  for  w  6  W 
may  be  understood  more  clearly  as  follows. 

(1)  Find  wi ,  the  longest  generator  in  G  which  is  also  a 
left  segment  (a  segment  mapping  from  the  same  initial  time 
as  w)  of  w. 

(2)  Repeat  for  what  remains  of  w  after  the  left  segment  wi 
has  been  removed.  If  the  process  stops  after  n  repetitions, 
then  w  =  w1»w2»...wn  and  {wi , . . .wn}  is  the  mis 
decompos i t i on . 

It  is  not  true  that  mis  decompositions  always  exist 
(that  is,  that  the  process  outlined  terminates).  Zeigler 
presents  a  theorem  for  conditions  guaranteeing  mis 
decompos i t i ons . 

If  G  generates  W  and  for  each  w  6  W  an  mis 
decomposition  by  G  exists,  G  is  said  to  be  an  admissible  set 
of  generators  for  W. 

The  above  discussion  of  operations  on  segments 
accomplishes  the  ability  to  move,  combine,  and  decompose 
them  canonically  under  stated  conditions.  This  framework  may 
be  used  in  describing  the  hierarchy  of  system 
speci f i cat i ons . 
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2.4  Hierarchy  of  System  Specifications 

There  is  a  heirarchy  of  specifications  for  systems  of 
increasing  power.  The  more  is  Known  about  a  system,  the  more 
detailed  and  definite  is  the  description  that  can  be 
specified.  At  the  level  where  a  system  is  Known  only  by  its 
behavior,  an  I/O  relation  observation  can  be  formulated. 

2.4.1  I/O  Relation  Observation 

At  the  most  basic  experimental  level,  a  system  is  Known 
by  its  input/output  pairs.  An  I/O  relation  observation  is  a 
structure  <T,X,W,Y,R>  where 
T  is  a  time  base 
X  is  an  input  value  set 
Y  is  an  output  value  set 

W  is  an  input  segment  set,  a  subset  of  (X,T) 

R  is  a  relation,  a  subset  of  W  x  (Y,T) 

where  (w,p)  6  R  =>  domain(w)  =  domain(p).  (that  is,  w 
and  p  are  input  and  output  segments  over  the  same 
observation  interval). 

R  is  the  collection  of  I/O  pairs,  and  the  structure  is 
a  purely  behavioral  description.  Nothing  is  said  about 
internal  state  of  the  system  or  about  means  of  generating 
output . 
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2.4.2  I/O  Function  Observation 

The  set  R  may  be  partitioned  if  different  initial 
states  of  the  system  are  related  to  different  sets  of  I/O 
pairs.  An  I/O  function  observation  is  a  structure 
<T,X,W,Y,F>  where  T,X,W  and  Y  are  as  before  and  F  is  a  set 
such  that 

f  6  F  ==>  f  is  a  subset  of  Wx(Y,T)  and  f  is  a 
function.  (If  p  =  f(w),  then  domain (p)  =  domain(w)). 

Each  f  refers  to  a  distinct  initial  state.  The 
collection  of  all  I/O  pairs  generated  by  all  f '  s  is  the  set 
R.  An  I/O  relation  observation  may  be  obtained  from  an  I/O 
function  observation,  then,  but  the  reverse  is  not  true  as 
information  about  initial  states  is  needed  for  an  I/O 
function  observation. 

2.4.3  Input  Output  System 

An  I/O  system  is  a  structure 
S  =  <T,X,W,Q,Y , f ,g> 

where 

T  is  a  time  base 
X  is  the  input  value  set 

W  is  a  subset  of  (X,T),  the  input  segment  set. 

Q  is  the  state  set 

Y  is  the  output  value  set 

f  is  the  state  transition  function 

g  is  the  output  function 

subject  to  the  constraints 


X 
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(1)  W  is  closed  under  composition,  that  is,  for  all 

contiguous  segments  wl ,w2  6  W,  w1*w2  6  W. 

(2)  f  and  g  are  deterministic  mappings 

f:  Q  x  W  -  -  >  Q 
g:  Q  -->  Y 

(3)  f  obeys  the  composition  property,  that  is,  for  every 
contiguous  pair  of  segments  wl ,w2  6  W 

f(q,w1»w2)  =  f ( f (q,w1 ) ,w2) 

The  significance  of  constraint  (1)  is  the  necessity  for 
compound  I/O  experiments  to  be  allowed  in  a  system. 
Constraint  (2)  embodies  the  requirement  of  a  deterministic 
system  that,  given  an  initial  state  and  a  particular  input 
segment,  the  system  will  inevitably  undergo  transition  to 
the  same  final  state.  Constraint  (3)  states  that  the 
fragmentation  of  an  input  segment  does  not  affect  the  final 
state  reached  through  operation  of  the  transition  function 
throughout  the  segment.  The  constraints  formalize  the  ideas 
intuitively  necessary  for  a  deterministic  dynamic  system, 
those  of  re-initialization,  repetition,  interruption  and  re¬ 
starting  of  an  experiment. 

Given  the  structure  S  it  is  possible  to  generate  the 
behavior  of  a  system.  That  is,  it  is  possible  to  obtain  an 
I/O  function  observation  and  an  I/O  relation  observation 
which  are  complete  in  the  sense  of  containing  all  possibly 
observable  I/O  pairs. 
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2.4.4  Iterative  Specification 

An  iterative  specification  of  a  system  is  a  structure 
G  =  <T , X , W , Q , Y , f , g> 

where  T,  X,  Q,  Y,  and  g  are  as  before,  and  W  is  a  subset  of 
(X,T),  the  input  segment  generator  set,  and  f:  Q  x  W  -->  Q 
is  the  single  segment  transition  function,  subject  to 

(1)  W  is  an  admissible  set  of  generators. 

(2)  f  has  the  composition  property. 

The  iterative  specification  structure  gives  only  the 
set  generating  the  input  segment  set  for  the  cor respondi ng 
system,  and  the  transition  function  for  each  segment  in  the 
generating  set.  Thus  if  w  has  an  mis  decomposition  by  W  of 
wl »w2»w3» . . . *wn ,  then  the  state  of  the  system  after  w  is 
given  by  applying  the  single  segment  transition  function  to 
each  wi  in  turn,  i.e.  if  initial  state  is  q(0),  then  let 

q ( 1 )  =  f(q(0),w1) 
q ( 2 )  =  f (q ( 1 ) , w2 )  etc . 

and  the  final  state  reached  is  q(n)  =  f(q(n-1),wn) 

The  utility  of  an  iterative  specification  is  twofold. 
First,  it  is  a  minimal  description  since  the  single  segment 
transition  function  operates  on  segments  which  may  occur  in 
the  decomposition  of  many  input  segments.  Second,  a 
transition  function  operating  iteratively  to  achieve  a 
compound  effect  is  useful  in  bridging  the  gap  between  model 
description  and  simulation  program. 

A  t ime- i nvar i ant  I/O  system  specification  can  be 
obtained  from  an  iterative  specification  by  compounding  the 
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single  segment  transition  function.  An  iterative 
specification  contains  more  information  in  that  it  is  Known 
that  the  input  segment  set  can  be  decomposed  and  the 
transition  function  can  be  represented  in  terms  of  action  on 
single  segments. 

2.4.5  Structural  Descriptions 

Zeigler  gives  a  means  of  describing  systems  which  can 
be  specified  as  interacting  networks  of  component  sub¬ 
systems.  The  descriptions  consist  of  specifications  of  the 
component  sub-systems  together  with  an  expression  labelling 
the  components  and,  for  each  one,  showing  those  others  which 
it  affects,  that  is  whose  input  segments  it  generates.  Such 
structural  specifications  correspond  most  closely  to 
informal  model  descriptions  in  terms  of  interacting 
components . 

Translation  of  a  specification  at  the  highest  level 
into  one  at  the  lowest  is  the  formal  equivalent  of  computer 
simulation  of  a  model.  Step-by-step  successive  states  and 
outputs  are  obtained  using  an  iterative  description,  and 
thus  behavioral  data  is  collected. 

2.5.  Model  Validity 

A  “model"  may  be  thought  of  as  one  system  which 
attempts  to  characterise  some  of  the  features  of  another 
system.  Successful  modelling  depends  on  the  establishment  of 
relationships  between  such  systems.  The  validity  of  a  base 
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model  in  representing  the  real  system,  of  a  lumped  model  in 
representing  the  base  model,  and  of  the  simulation  program 
in  implementing  the  lumped  model  must  all  be  investigated 
via  such  relationships. 

At  each  of  the  levels  of  specification  of  systems  a 
11  preservat  i  on  relation"  or  "morphism"  may  be  stated  which 
establishes  a  cor respondence  between  features  of  two  systems 
described  at  that  level.  To  the  hierarchy  of  specification 
levels  corresponds  a  hierarchy  of  morphisms.  This  hierarchy 
is  strict  in  that  lower  level  morphisms  do  not  imply  higher 
level  ones.  For  example,  systems  which  display  the  same  Kind 
of  input-output  behavior  may  generate  this  behavior  in 
entirely  different  manners. 

At  a  particular  level  of  specification,  morphisms  may 
be  described  between  systems  SI  and  S2 .  Let  SI  be  larger 
than  S2  in  the  sense  that  SI  is  capable  of  doing  anything 
that  S2  can,  but  not  vice  versa.  A  mapping  h  from  SI  to  S2 
is  called  a  "surjective"  or  "onto"  map,  and  ensures  that  all 
the  behavior  of  S2  is  being  covered  by  SI.  A  mapping  from  S2 
to  SI  \s  called  an  "injective"  map,  and  embeds  S2  in  SI,  or 
in  other  words,  isolates  the  appropriate  part  of  SI  to 
relate  to  S2 . 

2.5.1.  I/O  Relation  Observation  Morphism 

Let  SI  be  given  by  <T 1 , X 1 , W 1 , Y 1 , R 1 > .  Let  S2  be  given  by 
<T2 , X2 , W2 , Y2 , R2> .  Let  h:  W2  -->  W1 .  h  is  injective  and  can 
[30  referred  to  as  an  encoder.  If  w2  6  W2  is  applied  to  S2 , 
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then  h(w2)  =  wl  6  W1  tells  which  input  segment  to  apply  to 

51  . 

Let  K:  ( Y 1 , T 1 )  -->  (Y2,T2).  K  is  surjective  and  can  be 
referred  to  as  a  decoder.  An  output  segment  pi  6  ( Y 1  ,  T 1  )  is 
Known  to  represent  K ( p 1 )  =  p ( 2 )  6  ( Y 2  ,  T 2 ) . 

The  notion  that  the  I/O  behavior  R2  of  S2  can  be  found 
in  the  I/O  behavior  R1  of  S2  can  be  formalised  as  follows: 

An  "I/O  Relation  Observation  Morphism"  from  SI  to  S2  is 
a  pair  (h,K)  such  that 

1 .  h :  W2  - - >  W 1 

2.  K:  ( Y 1 , T 1 )  -->  (Y2,T2)  and  is  onto. 

3.  for  every  (w2,p2)  6  R2  there  exists  a  (w1,p1)  6  R1  such 

that  wl  =  h(w2)  and  K ( p 1 )  =  p2 . 

Thus  every  I/O  segment  pair  of  S2  has  its  corresponding 
pair  for  SI  and  the  mapping  of  one  pair  to  the  other  is 
established.  If  input  w2  is  applied  to  S2,  the  behavior  of 

52  may  be  studied  by  observing  the  behavior  of  SI  on 
application  of  h(w2),  and  decoding  the  result. 

2.5.2  I/O  Function  Observation  Morphism 

An  I/O  function  observation  morphism  from 
SI  =  <T 1 , X 1 , W 1 , Y 1 , F 1 >  to  S2  =  <T2 , X2 , W2 , Y2 , F2>  is  a  pair 
(h,k)  where 

1 .  h:  W2  -->  Wl 

2.  K:  ( Y 1  , T 1  )  -->  (Y2.T2) 

3.  For  each  f2  6  F2  there  is  an  fl  6  FI  such  that  f2  =  K*g; 

that  is,  for  all  w2  6  W 2,  f 2 ( w2 )  =  K ( f 1 ( h ( w2 ) ) ) . 
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The  effect  of  applying  w2  to  S2  may  be  studied  by 
coding  it  to  obtain  wl  =  h(w2),  applying  this  to  SI  and 
obtaining  pi  =  f 1 ( w 1 )  =  f1(h(w2))),  and  decoding  to 

p2  =v  K ( p 1 )  =  K( f 1 (h(w2 ) ) ) . 

2.5.3  I/O  System  Morphism 

A  system  morphism  from  a  system 

SI  =  <T 1 , X 1 , Wl , Q1 , Y 1 , f 1 ,g1 >  to  S2  =  <T2 , X2 , W2 , Q2 , Y2 , f 2 , g2> 
is  a  triple  (h,j,k)  where 

1 .  h:  Wl  -->  W2 

2.  j:  Q  -->  Q2  where  Q  is  a  subset  of  01,  and  j  is  onto. 

3.  k:  Q1  -->  Y2  and  is  onto. 

and  for  all  q  6  Q,  w2  6  W 2, 

4.  j ( f 1 ( q , h ( w2 ) ) )  =  f  2 ( j ( q ) , w  2 )  (transition  function 

preservat i on ) . 

5.  k ( g 1 ( q ) )  =  g2 ( j ( q ) )  (output  function  preservat ion ) . 

Under  a  system  morphism,  there  is  not  only  behavior 
preservation  but  also  there  are  structural  relationships 

between  the  state  spaces  and  the  means  by  which  state 

changes  and  behavior  are  generated. 

If  T 1  =  T2 ,  XI  =  X2 ,  Wl  =  W 2,  Y1  =  Y2  then  SI  and  S2 
are  called  "compatible".  Compatible  systems  SI  and  S2  are 
homomorphic  images  of  each  other  if  (h,j,k)  is  a  system 

morphism,  h  and  k  are  identity  maps,  and  Q  =  Q1 ;  j  is  said 

to  be  a  homomorphism.  If  j  is  one-to-one  as  well  as  onto  it 
is  an  isomorphism.  In  essence,  two  systems  are  homomorphic 
if  their  state  sets  map  completely  into  each  other  and  if 
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the  action  "encoding  (SI  to  S2)  followed  by  transition  (in 
S2 ) "  results  in  the  same  state  as  the  action  "transition  (in 
SI)  followed  by  encoding  ($1  to  S2) " , 

Systems  obeying  a  system  morphism  can  be  shown  to  obey 
behavioral  morphisms  as  well. 

2.5.4  Iterative  Specification  Morphism 

A  "specification  morphism"  is  a  triple  (h,j,K)  defined 
similarly  to  a  system  morphism,  with  the  differences  that  h 
maps  S2; s  generator  set  into  SI,  and  that  transition 
function  preservation  need  only  be  checked  for  the  generator 
set  of  the  smaller  system.  It  may  be  shown,  by  expanding 
each  iterative  specification  into  its  system  specifications 
through  compounding  the  single  segment  transition  functions, 
that  a  specification  morphism  implies  a  system  morphism. 
Thus  one  of  the  major  advantages  of  iterative  specifications 
is  the  ability  to  check  preservation  relationships  between 
them  without  expansion  to  system  form  and  exhaustive 
checking  of  state  transitions  for  all  members  of  the  input 
segment  set . 

2.5.5  Structured  Specification  Morphisms 

A  structured  specification  is  one  in  which  labels  and 
i nterconnect i ons  are  given  to  subsystems  of  a  large  system. 
Structural  morphisms  of  two  kinds  are  discussed  by  Ziegler 
(p.  276).  A  "weak  structure  morphism"  is  a  homomorphism 
between  the  state  structures  of  two  structured  systems,  but 
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while  the  global  transition  functions  of  the  systems  are 
related,  local  transition  functions  of  the  subsystems  need 
not  be.  (Typically  subsystems  of  the  larger  system  SI  are 
aggregated  into  blocks,  and  each  block  is  a  subsystem  of 
S2).  In  order  to  obtain  a  "strong  structure  morphism" 
between  systems,  it  is  necessary  to  ensure  that  local 
interaction  of  component  subsystems  is  preserved. 

The  morphisms  thus  constructed  play  a  fundamental  part 
both  in  consideration  of  model  (and  simulation)  validity  and 
in  facilitation  of  model  construction.  The  next  section 
presents  Ziegler's  twelve  fundamental  postulates  (Zeigler, 
1976)  for  successful  modelling  and  simulation,  and  these 
postulates  illustrate  the  utility  of  morphisms  in 
constructing  valid  models.  Ziegler  discusses  using  morphisms 
to  simplify  the  modelling  task  in  later  chapters  of  his 
book;  as  an  example,  large  systems  may  be  decomposed  into 
manageable  subsystems  provided  structural  morphisms  are 
mai ntai ned . 

High  level  morphisms  imply  lower  level  ones;  the 
significance  of  this  is  that  if  systems  obey  a  high  level 
morphism,  it  is  known  that  they  are  behavioral ly  related 
without  the  need  for  generating  behavior  and  comparing  it. 

2.6  The  Twelve  Postulates 

Zeigler  devotes  a  chapter  to  the  statement  and 
discussion  of  his  twelve  fundamental  postulates  (Zeigler, 
1976,  Chapter  11).  The  postulates  are  paraphrased  below. 
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Postulate  1:  There  exists  a  real  system  R  which  is 

identified  as  a  set  of  potentially  acquirable  data. 

Postulate  2:  There  exists  a  base  model  B  that  structurally 
characterizes  this  set  of  potentially  acquirable  data.  B  is 
a  t ime- i nvar i ant  transition  system  (system  without  output) 

B  =  <T,X, ( X , T ) ,  Q  ,  f  > 

Postulate  3:  There  exists  a  set  of  experimental  frames  { E } 
restricting  experimental  access  to  R. 

Postulate  4:  An  experimental  frame  E  is  a  structure 
E  =  <W,Y,g,V>.  That  is,  a  subset  W  of  (X,T)  of  input 
segments  applicable  in  E,  a  set  Y  of  output  values 
observable  in  E,  an  output  function  g  which  maps  base  model 
states  to  outputs,  and  a  set  V  of  valid  output  (determined 
by  experimental  conditions  of  E). 

Postulate  5:  The  real  system  observed  within  the 
experimental  frame  is  structurally  characterized  by  the  base 
model  in  the  frame,  denoted  B/E. 

Postulate  6:  B/E  =  <T , X , W , Q , Y , f , g>  where 
T  is  the  time  base  of  the  base  model 
X  is  the  input  value  set 

W  is  the  set  of  input  segments  allowed  under  E 
Q  is  the  set  of  base  model  states  which  lie  within  valid 
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range  for  E  together  with  0,  a  symbol  for  states  outside  E. 

Y  is  the  valid  output  set  together  with  0,  to  which  nonvalid 
outputs  are  collapsed. 

f:Q  x  W  -->  Q  is  the  transition  function  of  the  base  model 
where  it  results  in  valid  states  for  E;  all  other  states 
obtained  are  collapsed  to  0. 

g:Q  -->  Y  is  the  output  function  of  E  unless  applied  to  a 
state  not  valid  for  E,  in  which  case  it  maps  to  0. 

Postulate  7:  The  data  potentially  acquirable  by 
observations  of  the  real  system  within  experimental  frame  E 
are  identified  with  the  I/O  relation  of  B/E. 

Postulate  8:  The  real  system  R  is  the  set  of  data 
potentially  acquirable  by  observation  of  all  experimental 
frames . 

Postulate  9:  A  lumped  model  is  an  iterative  system 

specification  M. 

Postulate  10:  A  lumped  model  is  valid  for  the  real  system 
in  experimental  frame  E  if  its  1/0  relation  is  equivalent  to 
that  of  B/E. 

Postulate  11:  A  computer  is  an  iterative  system 
specification  C. 
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Postulate  12:  A  computer  C  is  a  valid  simulator  of  a  lumped 
model  M  if  there  is  a  specification  morphism  from  C  to  M . 

Postulates  10  and  12  define  what  is  meant  by  "valid 
model"  and  "valid  simulation",  and  the  limited  significance 
such  assessments  are  meant  to  have.  A  lumped  model  valid  for 
R  in  E  is  valid  in  terms  of  behavior;  its  own  Known 
structure  may  provide  insight  into  structural  properties  of 
the  base  model,  but  such  base  model  structure  cannot  be 
established  empirically.  However,  there  may  be  conditions 
under  which  structural  inference  about  B  may  be  made 
(Ziegler,  1976,  Chapter  15). 

2.7  Categories  of  Models 

Systems  may  be  categorized  according  to  whether  the 
time  base  is  discrete  or  continuous,  and  variable  sets  are 
finite,  discrete  or  continuous  (or  mixed). 

Models  in  common  use  most  often  fall  into  three 
categories;  di screte/di screte  (sequential  machines), 
continuous/discrete  ("discrete  event"  models)  and 
continuous/continuous  ("differential  equation"  models). 
Simulation  packages  for  models  in  the  latter  two  categories 
are  widely  available. 
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2.7.1  Discrete  Time  Models 

A  discrete  time  system  specification,  or  sequential 
machine,  is  a  structure 
M  =  <XM , QM , YM , f , g> , 

XM ,  QM ,  YM  subsets  of  the  integers.  Cor respondi ng  to  M  is  an 
iterative  specification  G  given  by 
G  =  <T , XM , W , QM , YM , f i ,gi> 

where 

T  =  I 

W  =  {w  |  w: (0, 1 )  -->  X} 

fi:  Q  x  W  -->  Q,  fi(q,w)  =  f ( q , w ( 0 ) ) 

gi  =  g 

Each  w:  (0,1)  -->  X  is  uniquely  defined  by  its  value 
w(0)  because  T  is  a  discrete  time  base.  Thus  the  generators 
in  W  associated  with  the  iterative  specification  are  in  a 
one-to-one  cor respondence  with  the  input  values  X.  Thus  the 
transition  of  the  system  under  f  in  response  to  a  sequence 
of  input  values  from  XM  is  equivalent  to  the  sequence  of 
transitions  fi  each  responding  to  a  single  input  from  XM. 

2.7.2  Discrete  Event  Models 

A  discrete  event  model  is  one  in  which  discrete  state 
changes  occur  at  arbitrary  time  intervals.  Events  are 
predictable  on  the  basis  of  previous  events  and  the  system's 
present  state,  and  therefore  can  be  scheduled. 

A  discrete  event  system  specification  is  a  structure 
M  =  <XM , $M , YM , f , g , t  > 
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where  XM,  YM  are  as  before,  and 
SM  is  the  set  of  sequential  states 
f  is  the  "quasi " -transi tion  function 
g  is  the  output  function 
t  is  the  time  advance  function  H 
such  that 

1)  t : s  -->  Rt  and  t(s)  is  the  time  the  system  is  to  remain 
in  state  s . 

2)  Let  QM  =  { (s,e)  |  s  6  SM ,  0<e<t(s)} 

Let  z  be  a  special  symbol  not  in  XM.  Then  f: 
QM  x  (XU  {z} )  -->  SM  such  that  f(s,e,z)  =  fz(s)  for  all 
(s,e)  6  QM ,  where  fz:SM  -->  SM . 

3)  g : QM  -->  YM 

Change  in  the  system  is  caused  either  by  completion  of 
time  t ( s )  and  arrival  of  the  next  (scheduled)  event,  or  by 
the  arrival  of  an  (unscheduled)  external  event.  If  no 
external  event  occurs  in  t(s),  symbolized  by  z,  then 
transition  is  to  the  next  state  fs(x).  If  an  external  event 
x  arrives  after  state  s  has  endured  for  time  e,  and  e<t(s), 
then  a  new  state  f(s,e,x)  is  immediately  initiated.  This  new 
state  depends  on  e,  and  thus  the  full  state  set  QM  includes 
e . 

An  iterative  form  of  the  system  specification  is 
constructed  by  considering  the  generating  segments  and 
associated  transitions.  Let  Wz  be  the  set  of  segments  of 
finite  length  where  no  external  event  occurs  and  W x  be  the 
set  of  segments  of  finite  nonzero  length  where  no  external 
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event  occurs  except  at  zero.  Then  the  generator  set  Wg  is 
the  union  of  Wz  and  Wx .  The  transition  function  can  then  be 
*  constructed  by  considering  actions  on  the  generating 
segments.  Inclusion  of  the  variable  e  in  the  state  set,  to 
show  elapsed  time  in  state  s,  means  that  in  the  iterative 
form  the  time  advance  function  can  be  dropped  and  a 
continuous  time  base  established.  Simulation  of  such  an 

iterative  form  is  then  straightforward . 

2.7.3  Differential  Equation  System  Specification 

Many  models  of  natural  phenomena  are  specified  in  terms 
of  instantaneous  rates  of  change  in  the  state  variables: 
that  is,  state  changes  are  not  prescribed  directly  but  must 
be  computed  as  satisfying  differential  equations. 

A  differential  equation  specified  system  is  a  structure 
D  =  <X,Q,Y,fp,g> 

where 

X  is  the  input  value  set 
Q  is  the  state  set 
Y  is  the  output  value  set 
fp  is  the  rate  of  change  function 
g  is  the  output  function 
subject  to 

(a)  X,  Q  and  Y  are  real  vector  spaces 

(b)  deterministic  responses 

fp:  Q  x  X  •■>  Q 
g:  Q  -->  Y 
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In  order  to  convert  this  specification  into  iterative 
form  it  is  necessary  to  assume  a  continuous  time  base  T,  and 
as  input  segment  generators,  the  bounded  continuous 
functions  from  finite  intervals  of  T  into  X.  Then  fp  can  be 
interpreted  in  a  manner  specifying  the  single  segment 
transition  function. 

Let  w:  <0,t1>  -->  X  and  q  6  Q  be  a  given  input  segment 
and  a  state.  Then  a  segment  P:  <0,t1>  -->  Q  is  a  solution 
associated  with  w  and  q  if  P ( 0 )  =  q  and 
dP(t)/dt  =  f p ( P ( t ) , w ( t ) )  for  each  t  6  <0,t1>.  If  exactly  one 
solution  exists  for  each  pair  (w,g),  the  specification  is 
that  of  a  deterministic  system  whose  state  trajectories  are 
given  by  the  solutions.  (Conditions  on  fp  which  ensure  the 
existence  and  uniqueness  of  solutions  are  Known). 

The  transition  function  in  the  iterative  specification 
of  a  unique  solution  differential  equation  system  is  given 
by  the  mapping  carried  out  by  the  solution  taking  a  given 
initial  state  into  the  cor respondi ng  final  state  at  the  end 
of  the  observation  interval. 

2.8  Simulation  Techniques  and  Facilities 

The  modeller  must  take  a  model  specification  in 
iterative  form  and  translate  it  into  a  program  i nterpretab 1 e 
by  the  computer.  If  the  model  falls  into  the  discrete  event 
or  differential  equation  category,  the  task  may  be 
simplified  by  the  availability  of  simulation  packages.  The 
modules  of  such  packages  may  have  specification  morphisms 
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Known  to  connect  them  to  the  model  elements  they  represent. 
If  the  model  does  not  fall  into  one  of  these  categories,  the 
programming  task  is  increased  unless  the  modeller  happens  to 
have  access  to  locally  available,  previously  written 
programs  applicable  to  his  work.  Programming  in  a  general 
purpose  language  may,  however,  result  in  a  more  intimate 
control  over  the  process  of  simulation. 


2.8.1  Continuous  Time  Simulation 

Implementation  of  continuous  models  on  a  digital 
computer  necessitates  discretization  of  the  continuous 
variables.  In  the  case  of  differential  equation  models, 
discretization  of  the  time  base  implies  approximation  of 
differential  equations  by  difference  equations.  The 
experience  of  the  numerical  analyst  in  choosing 
approximation  techniques  for  a  given  differential  equation 
may  then  be  applied  to  maintaining  a  specification  morphism 
between  the  transition  function  in  continuous  time,  as  given 
by  the  solution  of  a  differential  equation,  and  the 
transition  function  in  discrete  time,  as  given  by  the 
solution  of  the  approximating  difference  equation  by 
numerical  integration  techniques. 

Numerical  integration  packages,  which  contain 
programmed  versions  of  integration  methods,  and  associated 
documentation  assisting  the  modeller  in  choosing  a  method 


suited  to  his  needs,  are  widely  available.  Techniques  for 
estimating  errors  of  approximation  are  known  for  such 
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methods.  Furthermore,  simulation  languages  for  continuous 
modelling  such  as  CSMP  and  DYNAMO  are  ava i 1 ab 1 e  wh i ch 
perform  the  continuous  time  to  discrete  time  conversion  for 
the  modeller.  Successful  use  of  such  languages  depends  on 
the  appropr i ateness  of  the  approximations  chosen  for  the 
particular  problem;  the  modeller  may  need  closer  control. 

2.8.2  Discrete  Event  Simulation 

A  discrete  event  strategy  which  translates  a  model  into 
a  program  is  applicable  to  continuous  or  discrete  time 
models  in  which  state  changes  occur  at  variable  but 
predictable  time  intervals.  Fundamental  to  discrete  event 
techniques  is  the  formation  and  management  of  a  "next  event 
list"  which  orders  tasks  with  respect  to  time  and  allows 
discretization  of  the  time  base  into  the  appropriate  time 
steps.  Such  models  assume  that  nothing  happens  to  component 
state  variables  except  at  these  time  steps,  so  that 
discretization  error  between  model  and  computer 
respresentat i on  is  not  a  factor.  The  problem  in  maintaining 
morphisms  between  model  and  program  in  the  discrete  event 
strategy  is  that  of  ordering  and  executing  the  "next  event 
list"  in  a  manner  which  preserves  the  interdependence  of 
model  components. 

Three  basic  strategies  for  discrete  event  programming 
are  in  common  use  (Fishman,  1973).  For  each  strategy 
packages  are  commercially  available. 

The  "event  scheduling"  strategy  involves  recognizing  an 
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event  as  a  change  in  the  state  of  any  model  component.  The 
“next  event  list'1  is  initialized  with  an  event  which  is 
carried  out  (the  state  change  implemented),  and  any 
consequent  changes  in  state  of  any  component  entered  in  the 
appropriate  place  on  the  list.  Conflicts  arising  from  events 
scheduled  at  the  same  time  must  be  resolved  with  regard  to 
dependence  of  model  components  for  the  particular  model.  The 
simulation  languages  which  implement  this  technique  are  GASP 
and  SUBSCRIPT  which  are  both  FORTRAN  based.  They  may  be  used 
when  a  model  suggests  the  event  scheduling  strategy. 

The  “activity  scanning"  strategy  involves  recognizing 
as  entities  all  activities  in  the  model  which  can  affect  the 
states  of  components,  and  associating  with  each  a  set  of 
conditions  (on  model  component  states)  under  which 
initiation  or  termination  of  the  activity  occur.  For  each 
activity  there  is  a  clock,  or  sequence  of  time  steps  at 
which  state  changes  occur.  Processing  proceeds  by  checking 
the  clocks  of  all  activities  and  advancing  time  to  the  most 
imminent.  The  activity  scanning  technique  has  not  been 
widely  used,  possibly  owing  to  the  fact  that  a  simulation 
language  using  this  approach  is  not  widely  available. 

The  “process  interaction"  strategy  involves  recognition 
of  a  process  as  an  entity,  where  a  process  is  a  sequence  of 
events  changing  the  state  of  a  single  component.  Execution 
of  a  particular  process  may  have  to  be  halted  to  allow 
contingent  execution  of  other  processes  in  order  to  preserve 
the  interdependence  inherent  in  the  model.  In  this  category 
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there  are  two  widely  used  simulation  packages  available, 
GPSS  (written  in  IBM/360  assembly  language)  and  SIMULA 
(writ ten  i n  ALGOL ) . 

SIMULA  is  an  ALGOL-based  simulation  language  which 
allows  the  conceptual  operation  in  parallel  of  blocks  of 
ALGOL  programs.  As  well  as  controlling  execution  of  these 
blocks  by  using  the  process  interaction  technique,  it  allows 
description  and  generation  of  processes  during  simulation. 
The  modeller  describes  an  activity  procedure  for  each 
component  of  the  model,  and  SIMULA  controls  interaction  and 
contingent  processes  in  the  resulting  simulation.  SIMULA  is 
a  powerful  tool  where  applicable. 

2 . 9  Summmary 

In  this  chapter  an  outline  of  the  theory  and  conceptual 
elements  in  the  process  of  constructing  models  and 
simulating  them  has  been  presented.  Models  and  their 
computer  implementations  can  be  presented  in  a  formal 
language.  Within  this  framework,  the  success  with  which  one 
system  is  represented  by  another  can  be  formally  assessed. 
Models  can  moreover  be  viewed  in  the  perspective  of  classes 
of  models  and  facilities  available  for  conversion  to 
iterative  computer  form. 

The  remainder  of  the  thesis  will  be  presented  in  a 
format  cor respondi ng  to  the  five  elements,  the  real  system, 
the  experimental  frame,  the  base  model  (or  rather  models), 
lumped  models,  and  computer  implementations.  Questions  of 
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modelling  validity  will  then  be  considered  in  the  context  of 
simulation  experiments. 
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CHAPTER  THREE 


The  Real  System 


Since  the  discovery  of  the  genetic  code  and 

understanding  of  the  mechanisms  of  its  translation  into 
chemical  activity  in  the  cell,  a  major  area  of  study  for 
geneticists  and  developmental  biologists  has  been  that  of 
differentiation  and  development.  Scientists  have  been 
seeking  to  discover  how  cells  of  common  ancestry,  with  the 
same  information  stored  in  their  DNA,  proceed  to  interpret 
that  information  differentially,  and  thus  develop  into 
distinct  cell  types  in  the  cooperatively  structured 

organism.  The  problem  as  presently  conceived  is  that  of 
discovering  how  a  cell  has  in  active  i nterpretat i on  at  a 
given  time  only  those  portions  of  its  DNA  relevant  to  its 
present  function,  and  how  its  present  situation  leads  to 
other  portions  of  the  DNA  being  active  in  later  progeny 
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Theor i es 

of 

development 

have 

both 

guided 

exper imentat ion 

and 

been  illuminated  by 

i  t . 

It  is 

inevitable  therefore  to  discuss  such 

theor i es 

(models)  at 

the  same 

t  i  me 

as 

describing  the  biological 

systems 

themse 1 ves . 

Thus 

the 

background  of 

biology 

and 

existing 

theory,  out 

of 

which  the  author' 

s  own  mode  1 s 

arise,  is 

presented  as 

the 

"  rea  1 

system" . 

3.1  Genetic  Information  and  Development 

The  programmed  development  of  the  multicellular 
organism  seems  to  be  the  result  of  an  extremely  complex, 
multilevel  control  structure.  Various  experimental 
techniques  attack  particular  levels  in  the  hierarchy  of 
control . 

At  the  single-cell  level,  bacteria  have  been 
extensively  studied,  and  chemical  mechanisms  of  control  over 
their  relatively  minute  store  of  DNA  have  been  identified. 
Portions  of  the  DNA  code  for  proteins,  while  other  portions 
are  sites  involved  in  controlling  production.  In  the  well- 
known  '  operon'  model,  in  order  for  certain  code  to  be  read 
for  active  production,  a  portion  of  code  before  it,  known  as 
the  operator,  must  be  in  a  suitable  state.  By  binding  to 
the  operator,  certain  molecules  obstruct,  or  in  other  cases 
promote,  reading  of  following  code.  In  addition,  reaction 
with  some  other  substance  may  prevent  binding  by  these 
molecules.  Molecules  which  affect  control  areas  of  the  DNA 
may  either  be  extrinsic  to  the  bacterium  (thus  its 
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metabolism  can  adjust  to  its  chemical  environment),  or 
intrinsic,  that  is,  they  themselves  are  coded  for  in  the 
DNA.  In  this  way  the  production  of  one  metabolite  can 
affect  production  of  another,  or  a  critical  concentration  of 
a  metabolite  can  affect  its  own  production.  In  the 
bacterium,  therefore,  a  picture  emerges  of  control  at  the 
chemical  level,  with  several  types  of  feedback. 

In  multicellular  organisms  DNA  organization  is  thought 
to  be  immensely  more  complicated.  The  quantity  of  DNA  is 
orders  of  magnitude  larger.  Cellular  activities  are 
correspondingly  more  complex,  and  the  identification  of 
specific  chemical  mechanisms  in  control  becomes  very 
difficult.  Geneticists  using  mutation  experiments,  however, 
have  demonstrated  that  only  a  fraction  of  the  DNA  can  be 
directly  implicated  in  the  production  of  detectable 
substances  in  an  organism.  Furthermore,  some  mutations 
cannot  be  associated  directly  with  these  structural  regions 
of  the  DNA.  Since  such  mutations  frequently  have  serious 
consequences,  they  are  perhaps  the  result  of  changes  in 
control  portions  of  the  DNA  which  affect  future  development. 
It  has  also  been  demonstrated  that  those  portions  of  the  DNA 
which  appear  to  be  implicated  in  control  contain  reiterated 
code.  It  has  been  suggested  that  such  repetition  could  be  a 
means  of  timing  events  in  the  balanced,  se 1 f -regu 1  at i ng 
development  of  the  organism. 

A  model  showing  possible  levels  of  control  coded  for  in 
the  DNA  of  a  multicellular  organism  has  been  proposed  by 
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Britten  and  Davidson  (1969).  In  this  model,  batteries  of 
producer "  genes  are  activated  by  cascading  levels  of 
controlling  genes  which  in  turn  are  initiated  by  the  product 
of  a  single  "sensor"  gene.  Different  batteries  can  interact 
at  levels  below  the  sensor  gene,  if  products  of  genes  in  one 
battery  affect  the  activity  of  genes  in  another. 

The  analysis  of  DNA  structure  does  not  give  a  complete 
appreciation  of  the  mechanisms  of  development  and  cell 
differentiation,  however.  Development  is  above  all  a 
dynamic  process.  The  DNA  forms  the  data  for  a  program  which 
accesses  that  data  in  a  sequence  dependent  upon  the  results 
of  the  part  of  the  program  already  carried  out.  The  DNA  of 
a  cell  with  its  "active"  and  "inactive"  portions  may  be 
thought  of  as  the  cell's  "state".  Development  is  the  action 
of  a  "transition  function"  taking  a  cell  and  its  progeny 
from  one  state  to  another.  In  order  to  describe  the 
process,  one  of  two  methods  can  be  used.  In  the  first,  a 
"bottom-up"  method,  the  following  steps  would  be  taken: 

(a)  all  elements  in  the  state  of  a  cell  would  be  identified; 
that  is,  all  of  the  relevant  functional  products  of  the 
organism  would  have  to  be  listed; 

(b)  all  external  influences  on  (inputs  to)  a  cell  would  have 
to  be  i dent i f i ed ; 

(c)  a  "snapshot"  of  the  organism  for  each  state  occurring  in 
all  individual  cells  would  have  to  be  taken; 

(d)  with  all  of  the  above  data,  the  nature  of  the  transition 
function,  and  thus  of  the  developmental  process,  could 
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(possibly)  be  inferred. 

The  bottom-up  method  described  above  is  the  classic 
means  of  obtaining  an  operational,  dynamic  picture  of  a 
process  from  successive  static  descriptions.  It  is 
obviously  an  extremely  limited  approach,  both  practically 
and  conceptually,  in  the  case  of  complex  multicellular 
organisms.  For  that  reason  experimental  techniques  and 
control  models  of  a  "top-down"  nature  have  developed.  In 
this  second  approach,  the  following  steps  would  be  taken: 

(a)  a  developmental  phenomenon  would  be  identified  in  the 
organism  being  studied; 

(b)  any  regularly  occurring  events,  or  physical  patterns, 
which  might  be  associated  with  the  phenomenon  would  be 
observed ; 

(c)  experiments  designed  to  disturb  such  events  or  patterns 
would  be  performed  in  order  to  establish  possible  causal 
re  1  at i onshi ps ; 

(d)  once  such  causal  relationships  had  been  established, 
results  would  be  compared  to  others  obtained  in 
developmental  studies  to  test  their  generality  and 
plausibi 1 i ty . 

Two  instances  of  the  fruitful  use  of  a  top-down 
approach  will  be  given. 

In  a  well-known  embryo  1 og i ca 1  phenomenon  called 
i nduc t i on ,  when  a  certain  kind  of  tissue  (body  of  cells  with 
specialized  function)  grows  to  within  a  set  distance  of 
another  tissue  of  different  type,  a  change  is  observed  in 
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the  second  cell  type.  If  growth  to  the  critical  distance  is 
prevented,  no  change  occurs.  The  implication  is  that  some 
product  *  of  cells  in  the  first  tissue  triggers  a  chain  of 
events  in  the  second  type  which  results  in  a  change  in  the 
genetic  activation  of  the  second  type.  In  terms  of  an 
overall  picture  of  development,  such  phenomena  show  that  it 
is  important  to  think  in  terms  of  populations  of  cells  and 
interactions  between  them. 

As  a  second  example,  it  has  been  shown  that  the 
position  of  a  cell  in  a  Drosophila  imaginal  disc  (see 
Section  3.2)  determines  the  cell's  future.  The  implication 
is  that  cells  have  some  means  of  measuring  their  position  in 
a  population.  Again  the  significance  of  populations  of 
cells  and  signals  between  them  is  indicated. 

The  strengths  of  a  top-down  approach  are  that  it  does 
not  attempt  to  look  at  all  of  the  cell's  possible  activities 
at  once,  but  concentrates  on  a  certain  sequence  of  events, 
that  it  can  keep  as  many  cells  (as  much  of  the  organism)  as 
desired  in  view,  and  that  since  it  concentrates  on  change  it 
is  dynamic  in  approach. 

It  should  be  emphasized  that  the  top-down  type  of 
approach  advocated  here  is  a  "difference"  method.  It  does 
not  try  to  describe  a  cell's  state  of  activation  in  entirety 
but  rather  char acter i zes  differences  between  cells.  Studies 
of  induction  are  not  concerned  with  all  products  of  the 
cells  but  only  those  which  are  changed  by  the  inductive 
event.  The  characterization  of  populations  of  cells  in 
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terms  of  their  developmental  history,  then,  can  be  given  in 
terms  of  changes  in  behavior.  These  changes  and  differences 
between  cells  accumulate  as  development  progresses.  If  all 
differences  between  all  cells  at  different  stages  in 
development  were  to  be  catalogued,  together  with  the 
influences  causing  them,  a  complete  description  of  the 
genetic  apparatus  of  the  cell  would  be  achieved.  The  power 
of  a  difference  method  is  that  it  proceeds  from  simplicity 
(cells  which  are  initially  not  mutually  differentiated)  to 
complexity,  in  a  dynamic  manner,  without  trying  to  maintain 
a  universal  perspective. 

3.2  Disc  Studies  in  Drosophila 

In  the  larval  stage  of  development  of  the  fruit  fly 
Drosophila  small  buds  of  tissue  called  imaginal  discs  form. 
They  are  destined  to  become  the  external  limbs  and  organs  of 
the  adult  fly.  These  discs  can  be  manipulated 
experimentally,  cultured  in  the  abdomens  of  adult  flies, 
implanted  back  into  larvae  and  allowed  to  differentiate 
during  metamorphosis.  The  results  of  the  experimental 
manipulation  can  be  observed  in  the  abnormal  structures 
resulting  in  the  mature  host  fly.  From  the  extensive 
research  on  imaginal  discs,  a  few  pertinent  results  will  be 
described  here. 

3.2.1.  Transdetermination  and  Codes 

Hadorn  (1966)  has  shown  that  discs  undergo 
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determination  with  respect  to  their  future  organ  type  at  an 
early  stage,  long  before  visible  expression  of  this  type. 
All  progeny  cells  inherit  the  same  determination  as  their 
parent  cells.  Hadorn  found  that  with  repeated  culture  of 
discs  in  adult  abdomens,  a  small  fraction  of  these  discs 
change  their  state  of  determination  and  become  forerunner 
cells  for  another  disc  type.  This  phenomenon  is  Known  as 
transdetermi nation.  The  frequency  of  transdetermi nat ion 
between  the  various  disc  types  was  calculated  for  aabout 
6,000  implanted  fragments  observed  through  metamorphosis. 

Kauffman  (1973)  has  used  Hadorn' s  transdetermi nat ion 
data  as  the  basis  for  a  model  of  determination  in  terms  of  a 
number  of  bistable  switches.  The  idea  is  that 
transdetermi nat ion  occurs  because  the  disc-type  decision  is 
of  a  switch- like  nature,  and  decisions  can  be  changed  by 
resetting  switches.  Kauffman's  assignment  of  bistable 
switches  is  shown  in  Figure  3.1. 

The  frequency  of  transdetermi nat i on  should  be  inversely 
proportional  to  the  number  of  switches  which  are  in  a 
relatively  stable  state.  A  minimum  of  4  switches  is 
required  for  good  agreement  with  the  frequency  data.  The  1 
state  of  a  switch  is  assumed  more  stable  than  the  0;  thus 
frequencies  of  transdetermination  in  opposite  directions  are 
not  equal.  For  example,  Eye  to  Wing  happens  more  often  than 
Wing  to  Eye,  as  suggested  by  the  lengths  of  the 
cor respondi ng  arrows. 

Kauffman's  model  is  a  difference  model  of  the  type 
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FIGURE  3.1  Kauffman's  assignment  of  bistable  switches 
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discussed  in  the  previous  section.  It  is  not  a  picture  of 
activation  and  de- act i vat i on  in  the  DNA  cor respondi ng  to 
disc  type,  but  rather  a  minimal  description  of  differences 
among  cell  types  and  transformations  between  them.  The 
model  is  important  because  it  expresses  these  differences  in 
quantitative  terms.  Hamming  distance  between  the  codes  for 
different  discs  has  meaning  in  terms  of  genetic  similarity 
between  them. 

3.2.2.  Positional  Information  and  Gradients 

A  number  of  experiments  indicate  that  the  location  of  a 
cell  in  a  disc  has  more  significance  than  determining  the 
final  structure  of  the  organ.  For  example  Schubiger  (1971) 
has  found  differences  between  parts  of  the  disc  in  terms  of 
their  capacity  to  produce  duplicates  of  the  features  they 
determine,  when  fractured  from  the  rest  of  the  disc  and 
allowed  to  grow,  or  to  reproduce  features  which  were 
removed.  He  has  also  found  a  range  of  frequencies  of 
transdetermination  for  different  fragments.  So,  before  any 
question  of  final  differentiated  cell  type  arises,  cells 
have  some  way  of  sensing  where  they  are  in  the  disc;  that 
is,  the  cells  possess  positional  information. 

Bryant  (1971)  has  developed  the  idea  of  an  information 
gradient  in  cells.  In  situ  cuts  of  discs  into  two  pieces 
result  in  the  regeneration  of  the  entire  disc  from  one 
portion  and  mirror  image  formation  by  the  other,  judged  by 
the  regenerative  or  duplicative  effects  seen  in  the  adult 
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organism.  This  phenomenon  can  be  explained  in  terms  of  a 
gradient  of  regenerative  potential,  as  shown  in  Figure  3.2. 
The  portion  of  disc  containing  the  regenerative  hotspot  A 
will  regenerate,  while  the  other  portion  can  only  produce 
features  lower  than  its  own  highest  point  on  the  gradient. 
Region  A  has  been  located  on  the  disc. 

A  large  amount  of  evidence  points  to  this  capability  of 
cells  to  read  information  of  a  positional  nature.  As  one 
example,  if  a  portion  of  the  developing  leg  of  a  cockroach 
is  removed,  compensatory  growth  occurs.  This  can  be 
explained  in  terms  of  the  existance  of  a  gradient  down  the 
leg;  when  a  portion  is  removed,  a  discontinuity  is  caused  in 
the  gradient  which  acts  as  a  signal  to  cells  to  divide  in 
that  region. 

French,  Bryant  and  Bryant  (1976)  have  developed  a 
formal  model  for  pattern  regulation,  using  polar  coordinates 
to  describe  a  gradient  field  in  two  dimensions.  There  is 
direct  evidence  for  two  simple  rules  governing  regeneration 
in  cockroach  legs,  imaginal  discs  and  amphibian  limbs.  The 
"shortest  intercalation  rule"  stipulates  that  when  normally 
nonadjacent  positions,  either  in  the  radial  or  the  circular 
sense,  are  confronted  in  a  graft  combination  or  in  wound 
healing,  growth  occurs  at  the  junction  until  all 

intermediate  positions  have  been  intercalated.  (If  a  circle 
of  positions  is  interrupted,  missing  positions  on  the 
shortest  arc  joining  the  confronting  cells  are 

intercalated.)  The  "complete  circle  rule"  stipulates  that  an 
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FIGURE  3.2  Bryant's  gradient  model 
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entire  circle  of  positional  values,  at  a  certain  radial 
position,  may  undergo  distal  t r ansformat i on  to  produce  cells 
with  all  the  more  central  (distal)  positional  values.  That 
is,  the  circle  will  produce  all  values  contained  in  a  disc 
extending  to  that  radial  position.  The  authors  suggest  that 
since  the  model  can  explain  phenomena  in  diverse  organisms, 
it  may  be  generally  applicable  to  pattern  formation. 

Several  questions  are  open  concerning  gradients.  Some 
of  these  are: 

(a)  What  is  the  physical  basis  for  these  gradients?  Are 
they  as  simple  as  concentration  gradients? 

(b)  How  many  different  fields  of  positional  information  are 
needed  to  give  a  cell  its  future  developmental  instructions? 

(c)  What  shapes  of  fields  exist?  Are  they  rectangular, 
spherical,  or  perhaps  spiral  chemical  waves? 


3.2.3.  Compar tmenta 1 i zat i on 

Work  of  great  interest  has  been  done  in  the  last  few 
years  in  the  laboratory  of  Garci a-Bel 1 ido  and  his  co-workers 
(1973).  An  excellent  summary  of  this  work  can  be  found  in  a 
paper  by  Crick  and  Lawrence  ( 1975) ,  who  indicate  the 
importance  of  the  results  to  developmental  questions. 

A  mutation  in  a  particular  cell  is  called  'gratuitous' 
if  it  does  not  interfere  with  the  normal  development  of  the 
organism.  If  a  gratuitous  mutation  results  in  progeny  cells 
which  are  identifiable  against  a  background  of  non-mutant 
cells,  on  the  epidermis  for  example,  a  technique  known  as 
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clonal  analysis  can  be  used.  In  this  technique  lineages  of 
cells  are  followed  by  studying  the  descendants  of  a  number 
of  marked  cells  in  different  specimens.  Garci a- Be  1 1 i do  et 
al.  marked  cells  in  this  way  by  using  irradiation  of  the 
wing  disc  in  a  special  genetic  strain  of  Drosophila. 

The  earlier  a  cell  is  marked  in  the  developing  disc, 
the  larger  the  number  of  marked  descendants  there  will  be  in 
the  observable  mature  wing.  The  ratio  of  marked  to  unmarked 
cells  indicates  how  many  cells  there  were  at  the  time'  of 
irradiation.  If  there  were  n  cells  at  irradiation,  and  one 
is  marked,  then  1/n  of  all  cells  in  the  mature  wing  will  be 
marked.  The  patches  of  marked  cells  were  initiated  at 
various  times  in  the  development  of  the  disc.  They  show 
coherence,  and  the  startling  phenomenon  which  has  come  to  be 
known  as  compar tmenta 1 i zat i on .  Marking  cells  one  generation 
apart  leads  to  the  discovery  that  the  disc  undergoes  a 
number  of  subdivisions  into  distinct  areas  with  boundaries 
impervious  to  intrusion  by  progeny  cells  from  bordering 
areas.  A  cell  marked  just  before  such  a  division  may  have 
progeny  in  both  compartments  while  a  cell  marked  after  the 
division  does  not.  If  the  clone  of  marked  cells  on  the  wing 
encounters  a  boundary  the  patch  follows  it.  Thus  it  has 
been  discovered  that  the  lines  of  division  are  smooth  and 
invariant  from  specimen  to  specimen.  Furthermore,  the 
division  events  occur  in  a  set  sequence  and  are  inescapably 
implicated  in  the  developmental  program  of  the  organism. 

Investigations  in  other  discs, 
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have  since  shown  that  compar tmenta 1 ization  is  a  frequent 
event  in  development.  It  may  be  even  more  frequent  than  can 
be  seen  because  it  is  readily  recognisable  only  in  systems 
of  a  convenient  two-dimensional  nature  such  as  the  epidermis 
of  a  wing. 

Crick  and  Lawrence  see  three  possible  mechanisms  for 
the  establishment  of  compar tmenta 1  boundaries. 

1.  Partition  of  daughters:  each  cell  has  one  daughter  which 

ends  up  in  one  compartment  and  one  in  another. 

2.  Random  allocation:  cells  are  allocated  to  compartments 
at  random  according  to  some  statistical  measure  keyed 
to  the  size  of  the  resulting  compartments. 

3.  Geographical  partition:  founder  cells  of  one  compartment 

are  separated  from  those  of  another  by  the 

establishment  of  a  physical  line  of  division. 

The  first  two  mechanisms  seem  unlikely,  as  a  large 
amount  of  cell  movement  and  sorting  out  would  be  necessary 
to  allow  grouping  of  the  cells  of  a  compartment.  That  much 
movement  is  unlikely  in  view  of  the  coherent  clones  of 
marked  cells  which  have  been  observed.  The  third  mechanism 
appears  the  most  likely,  and  would  be  the  result  of  a 
phenomenon  of  field- like  nature.  Gradients  of  information 
across  the  disc  again  appear  necessary. 

Crick  and  Lawrence  propose  a  list  of  questions 
stimulated  by  compar tmenta 1 ization  studies.  They  are 
summarized  below. 

(a)  Are  all  divisions  binary?  If  they  are,  does  this  imply 
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a  binary  switching  mechanism  in  the  DNA  itself? 

(b)  Why  are  compartment  boundaries  smooth,  even  to  within  a 

cell  diameter?  Do  boundaries  lie  along  smooth  contours 
of  fields  of  information,  or  is  there  cell  sorting 
between  cells  of  distinct  characters,  or  a  combination 
of  the  two  mechanisms? 

(c)  When  do  subdivisions  cease?  Does  this  Kind  of 

phenomenon  continue  until  final  cell  differentiation 
itself,  or  are  there  other  mechanisms  at  work? 

(d)  How  are  cells  in  different  compartments  actually 

different?  Are  there  immediately  expressed  differences 
as  well  as  final  differentiated  cell  types  and  organ 
parts? 

(e)  What  is  the  connection  (if  any)  or  conflict  (if  any) 

between  compartments  and  gradients  of  information? 

(f)  Are  compar tmenta 1  decisions  i r revers i b 1 e?  How  stable 

are  genetic  decisions? 

(g)  What  is  the  significance  of  the  coincidence  of 

compar tmenta 1  boundaries  with  regions  of  uniformly 
transformat  ion  in  homeotic  mutants  (mutants  in  which 
one  tissue  is  uniformly  transformed  into  another  in  a 
specific  region)  (Kauffman,  1978)?  Can  the  homeotic 
genes  be  identified  with  Britten  and  Davidson's  (1969) 
sensor  genes? 

In  a  subsequent  paper  Lawrence  (1975)  has  attempted  to 
discover  the  degree  of  discontinuity  across  compartment 
boundaries.  The  cells  appear  to  be  functionally  distinct  in 
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some  but  not  all  respects.  Communication  of  a  chemical 
nature  across  boundaries  is  possible  but  may  be  selective. 

In  summary,  compar tmenta 1 i zat i on  appears  to  be  an 
essential  part  of  development  in  many  organisms.  It  appears 
to  be  the  basis  for  organization  of  populations  of  cells 
into  functionally  distinct  groups  and  populations  of 
differentially  communicating  cells.  The  connection  between 
compartments  and  gradient  fields  of  information  is  an 
important  focus  of  attention  in  present  studies.  One  of  the 
major  difficulties,  of  course,  is  that  the  nature  of 
gradients  remains  largely  obscure. 

Three  fundamental  ideas  have  been  presented  in  this 
discussion  of  issues  in  genetics.  They  are  char acter i zat i on 
of  the  problem  in  terms  of  difference  models  of  a 
quantitative  nature,  the  importance  of  positional 
information  in  development,  and  the  phenomenon  of 
compar  tmenta 1 i zat i on . 


3.3  The  Mechanics  of  Pattern  Formation 

A  fundamental  process  in  development  is  the  formation 
of  spatial  patterns  of  differentiated  tissues  from  initially 
homogeneous  tissue.  Individual  cells  in  a  specific  tissue 
have  differentiated,  but  formation  of  the  spatial  structure 
of  the  tissue  necessitates  a  super -ce 1 1 u 1 ar  mechanism.  It 
has  long  been  thought  (Child,  1941;  Waddington,  1962)  that 
"primary  fields"  of  morphogen  concentrations,  or  of  other 
spatially  varying  physical  parameters,  must  be  established 
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as  pre-patterns  for  tissue  structures.  The  genetic 
constitution  and  history  of  a  cell  provides  specificity  of 
response  to  positional  signals;  differential  response  of 
cells  of  similar  genetic  constitution  depends  on  positional 
specificity  extrinsic  to  the  genome  ( Garci a-Bell ido,  1975). 
Thus  topological  organisation  precedes  and  conditions 
specific  gene  activation. 

Investigations  of  various  i nterce 1 1 u 1 ar  mechanisms 
capable  of  establishing  spatial  patterns  have  been 
undertaken.  These  investigations  are  typical  of  the  "top- 
down"  approach  of  the  developmental  geneticist  in  that 
specific  interaction  of  morphogenetic  agents  and  the  DNA  is 
in  many  cases  unknown.  However  if  spatio-temporal  patterns 
of  such  agents  are  found  to  coincide  with  developmental 
structures,  a  causal  (or  at  least  influential)  relationship 
may  be  inferred. 

A  morphogenetic  field  may  be  defined  as  an  assembly  of 
functionally  coupled  cells  whose  development  is  under  the 
control  of  a  common  regulatory  process.  Establishment  of  a 
pre-pattern  in  such  a  field  implies  a  mechanism  by  which  an 
initially  homogeneous  distribution  of  a  morphogenetic  agent 
evolves  a  spatial  inhomogeneity  capable  of  initiating  and 
sustaining  a  pattern  of  differentiating  cells.  Several  such 
mechanisms  have  been  studied. 
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3.3.1  Active  Transport 

If  a  polarized  asymmetry  exists  in  each  of  the  cells  of 
the  field,  a  morphogen  can  be  actively  propogated  in  a 
particular  direction  through  the  field,  resulting  in  a 
gradient  of  concentration.  Boundary  cells  in  such  systems 
have  special  properties  sustaining  the  gradient  (Lawrence, 
1966;  Wi Iby  and  Webster,  1970).  A  similar  system  involving 
boundary  "sources"  and  "sinks"  and  unpolarized  cells,  has 
been  developed  by  Babloyantz  and  Hiernaux  (1974). 
Inhomogeneity  in  such  models  is  a  result  of  a  mechanism 
which  actively  opposes  the  natural  diffusion  of  molecules 
down  their  concentration  gradient  to  achieve  uniform 
concentration,  the  minimum  energy  configuration. 

3.3.2  Phase  Gradients 

If  the  cells  of  a  field  can  be  regarded  as  autonomous 
oscillators  which  emit  signals  of  short  duration  at  regular 
intervals,  and  the  signals  of  a  cell  are  propogated  by  other 
cells  after  a  delay,  a  wave  of  concentration  from  each  cell 
arises.  If  a  gradient  of  emission  intervals  is  pre-existent 
in  the  field  (implying  an  initial  polarity),  a  stable 
pattern  of  positional  information  can  result  (Goodwin  and 
Cohen,  1969).  Aggregation  in  cellular  slime  molds  is  an 
example  of  the  organisational  capability  of  oscillatory 
processes  (Nicolis  and  Prigogine,  1977). 
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3.3.3  Cellular  Contacts 

Positional  specificity  of  morphogen  concentration  can 
be  obtained  from  an  initially  polarized  field  without 
widespread  propogation  if  the  concentrat ion  level  in  a  cell 
is  affected  by  the  neighbour  in  the  polar  direction 
(Babloyantz,  1977)  or  neighbours  in  bipolar  directions 
( Garay ,  1 977 ) . 

3.3.4  Cell  Sor t i ng 

If  it  is  assumed  that  the  genetically  identical  cells 
of  a  field  are  in  some  sense  differentiated  with  respect  to 
some  adhesive  or  chemotaxic  property,  then  conceivably  a 
cell-sorting  mechanism  could  give  rise  to  the  assembly  of 
patterns  of  cells  by  autonomous  movements  of  cells  in  the 
field . 

All  of  the  above  mechanisms  involve  either  initial 
polarity  of  the  field  or  initial  asymmetry  of  cells,  or 
both,  and  thus  violate  strict  i nterpretat i on  of  the  initial 
homogeneity  of  the  field.  It  has  been  shown  that  processes 
naturally  occurring  in  fields  of  cells,  which  initially 
achieve  steady  state  (time  invariant),  spatially  uniform 
conditions,  may,  under  certain  influences  also  naturally 
present,  deviate  from  such  conditions  and  produce  a  non- 
uniform  stable  state.  An  initial  investigation  of  such 
situations  was  undertaken  by  Turing  (1952).  Recently 
several  theories  of  pattern  formation  based  on  similar 
phenomena  have  been  formulated  (Gmitro  and  Scriven,  1966; 
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Crick,  1970;  Gierer  and  Meinhart,  1972;  Meinhart,  1977; 
Kauffman  et  al.,  1978). 

3.3.5  Reaction-Diffusion  Kinetics 

Consider  a  field  in  which  local  conditions  are  defined 
by  the  concentrations  of  various  chemical  species.  Any 
i rregu 1 ar i t i es  in  these  concentrations  tend  to  be  damped  by 
diffusion;  that  is,  concentrations  tend  towards  spatial 
uniformity  by  transport  of  the  species  down  its 
concentration  gradient.  Diffusion  may  be  regarded  as  a 
consequence  of  the  second  law  of  thermodynamics,  which 
states  that  in  a  closed  system  entropy  increases 
monotonical ly  until  its  maximum  is  reached  at  "thermodynamic 
equilibrium".  Thermodynamic  equilibrium  in  the  field 
corresponds  to  steady  state,  spatially  uniform 
concent r at i on . 

I r regu 1 ar i t i es  in  concentrations  are  caused  either  by 
chemical  reaction  within  the  field  or  by  mass  exchange  from 
the  sur roundi ngs .  For  the  concentration  of  a  chemical 
species  X(r,t),  the  following  reaction-diffusion  equation 
can  be  established: 

d(x)/d(t)  =  F(X,Y,Z, . . . )+D( 1 )L(X)+M(x)  (1) 

That  is,  the  concentration  evolves  in  time  according  to: 

(a)  chemical  reaction,  F(X,Y,...)  (which  may  involve  other 
species  Y,Z,...);  (b)  diffusion,  which  may  be  approximated 

in  dilute  mixtures  by  an  independent,  constant  diffusion 
rate  D(1)  multiplying  the  Laplacian  operator  (sum  of  second 
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spatial  derivatives,  or  rate  of  flux  of  X  at  r);  (c)  mass 
influx  through  the  boundaries  M ( X ) . 

An  equation  of  this  type  is  classified  as  a  parabolic 
problem  in  that  the  time  derivative  is  of  first  order  and 
spatial  derivatives  are  second  order .  The  spatial  pattern 
of  concentration  of  X  will  depend  on  the  stability  of  the 
thermodynamic  equilibrium,  steady  state,  spatially  uniform 
solution  given  by  F ( XO , YO , ZO , . . . ) +M ( XO ) =0 .  The  term  M ( X ) 
may  be  considered  as  introducing  boundary  conditions;  M(x)=0 
on  the  boundary  corresponds  to  a  closed  system  with  no  flux 
across  the  boundary  (a  "Neumann"  condition),  while  M ( X )  of  a 
certain  value  on  the  boundary  (a  "Dirichlet"  condition) 
corresponds  to  an  open  system  (for  example,  "source"  and 
"sink"  problems).  Specified  boundary  conditions  for  the 
reaction-diffusion  equation  ensure  the  uniqueness  of  the 
steady  state  solution  XO. 

It  has  been  shown  (Turing,  1952;  Nicolis  and  Prigogine, 
1977)  that  for  certain  types  of  reaction  functions 
F(X,Y,Z,...)  the  thermodynamic  equilibrium  solution  XO  is 
not  stable,  and  per turbat ions  produce  stable,  spatially 
inhomogeneous  patterns  of  concentration.  Some  analytical 
techniques  described  below  may  be  used  to  examine  this 
s i tuat ion . 


1.  Bifurcation  Theory.  Bifurcation  theory  was 
Poincare  and  developed  by  other  mathematicians 
non-equilibrium  solutions  of  nonlinear 
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equations  (Nicolis  and  Prigogine,  1977). 

Suppose  in  equation  (1)  that  M(X)=0  and  F  is  governed 
by  a  set  of  parameters  p  (corresponding  to  rates  of 
reaction).  If  the  p  are  small  enough,  solution  of  the 
system  evolves  towards  thermodynamic  equilibrium  (XO)  as  t 
tends  to  infinity.  Now  for  critical  values  of  p  the  system 
can  evolve  to  a  new  steady-state  solution  due  to  instability 
of  the  thermodynamic  equilibrium.  Such  values  of  p  are 
called  "bifurcation  points".  The  purpose  of  the  theory  is 
to  develop  methods  to  (a)  locate  such  points  for  a  given 
system  and  (b)  construct  approximate  analytic,  convergent 
expressions  for  solutions  emerging  at  bifurcation  points. 
An  example  of  the  application  of  bifurcation  theory  to  a 
reaction-diffusion  system  may  be  found  in  Nicolis  and 
Prigogine  ( 1977)  . 

2.  Theory  of  Catastrophes.  In  further  development  of  the 
work  of  Poincare,  Thom  (1972)  has  developed  a  means  of 
classifying  systems  which  can  undergo  transitions  from  one 
steady  state  to  another.  Such  systems  must  be  expressible 
as  first  order  ordinary  differential  equations  with  time  as 
the  independent  variable.  The  diffusion  term  in  (1)  is  not 
present  explicitly;  all  spatial  dependence  is  expressed  by  a 
parametric  dependence  of  F  in  space  and  time.  The  time 
derivative  arises  from  a  "potential"  function 

dX/dt  =  -  dV ( X , Y , Z , . . . )/dX  (2) 

Steady  states  of  the  system  correspond  to  dV/dX=0  and  exist 
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if  V  does  not  depend  explicitly  on  t.  dV/dX=0  defines 
hyper sur faces  in  the  parametric  space.  Identification  of 
points  where  there  may  be  transition  between  steady  states 
becomes  a  topological  problem  of  identifying  regions  in 
which  a  curve  (particular  solution)  may  cross  from  one 
hypersurface  to  another.  An  example  of  the  application  of 
catastrophe  theory  to  a  one  species  reaction  system  may  be 
found  in  Nicolis  and  Prigogine  (1977). 

The  potential  importance  of  catastrophe  theory  in 
describing  phenonema  of  biological  order  has  been  emphasized 
by  Goodwin  (1973).  Severe  criticisms  of  some  attempted 
applications  have  been  made,  however  (Zahler  and  Sussman, 
1977);  it  appears  that  a  fundamental  problem  arises  out  of 
the  present  lack  of  quantitative  guidelines  to  its  practical 
application  in  problems  of  development. 

3.  Stability  Theory.  Behavior  of  a  reaction-diffusion 
system  such  as  (1)  under  perturbation  of  the  equilibrium 
solution  may  be  studied  by  investigation  of  the  stability  of 
that  equilibrium  configuration  XO.  Stability  can  be 
investigated  by  "Lyapunov's  Second  Method"  (Nicolis  and 
Prigogine,  1977),  in  which  a  functional  must  be  constructed 
of  the  concentrations.  The  functional  plays  the  role  of  a 
potential  in  the  vicinity  of  XO.  Asymptotic  stability  at  XO 
is  established  if  a  (local)  absolute  minimum  of  V  exists 
there.  Asymptotic  stability  is  the  property  that  any 
solution  in  a  small  neighbourhood  a  of  XO  at  time  t ( 0 ) 
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approaches  XO  as  time  tends  to  infinity. 

The  stability  of  XO  may  also  be  investigated  by  means 
of  linearization  of  the  system  in  the  neighbourhood  of  XO. 
This  method  has  been  investigated  in  detail  by  Gmitro  and 
Scriven  (1966),  is  discussed  by  Nicolis  and  Prigogine 
(1977),  and  has  been  practically  applied  to  specific 
problems  in  pattern  formation  by  Turing  (1952),  Meinhardt 
(1978),  Wi lby  and  Ede  (1976),  Kauffman  et  al.  (1978),  and 
others . 

The  steady  state  ( t i me- i nvar i ant ) ,  spatially 
homogeneous  (diffusion  term  zero)  solution  XO  of  (1)  obeys 

0  =  F ( XO , Y,Z, . .  . )+M(X0)  (3) 

Unsteady  states  can  be  represented  by  the  sum  of  steady 
state  concentration  XO  and  a  depature  from  the  steady  state, 
x;  that  is,  X=X0+x. 

Equation  (1)  can  the  be  written 

d/dt  ( XO  +  x ) =F ( XO  +  x , Y , Z ,  . .  . )+D(  1  )  L ( XO  +  x ) +M ( XO  +  x )  (4) 

Under  approppriate  smoothness  conditions,  F  and  M  can  be 
expanded  in  terms  of  a  Taylor's  series  in  x  around  XO; 
further,  if  x  is  small  enough,  terms  of  order  x**2  and 
higher  may  be  neglected  in  these  expansions.  Then  (4) 
becomes 

dXO/dt+dx/dt =F ( XO , y , z , . . . )+dF (X0)x/dx+D( 1  )  1  ( XO ) +D ( 1 ) 1 ( x )  + 

M ( XO ) +M ( x )  (5) 

By  (3)  and  dX0/dt=0  this  becomes 

dx/dt=dF (XO) x/dx+dM ( XO ) x+D ( 1 )L(x) 
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or  dx/dt=Kx+D ( 1 ) L ( x )  (6) 

This  is  a  linear  equation  which  may  be  solved  by  the 
methods  of  harmonic  analysis.  For  several  species 
X ( 1 ) , X ( 2 ) , X ( 3  )  ,  .  .  .  whose  reaction-diffusion  equations  are 
coupled,  K  is  the  Jacobian  matrix  of  F+M  with  respect  to 
x (  1  ) , x ( 2 ) , x ( 3 ) , . . .  evaluated  at  the  thermodynamic 
equilibrium  solution  of  the  system  of  equations. 

The  methods  of  harmonic  analysis  are  based  on  the 
fundamental  theorem  of  Fourier  which  states  that  any  spatial 
configuration  in  a  domain  may  be  expressed  as  a  weighted  sum 
of  elementary  shape  functions  appropriate  to  that  domain. 
The  weighting  factors  relate  the  time-dependent  pattern 
particular  to  x  to  the  elementary  patterns.  Thus  (x)  (the 
column  vector ( x ( 1 ), x ( 2 ),..., x ( n )) )  may  be  written 
( x ) = ( a ( K , n ) (A) (k,n) exp (L(k,n)t))F(k,r)  (7) 


where 

a)  F(k,r)  are  the  shape  functions  of  the  domain 

b)  L(k,n)  and  (A)(k,n)  are  the  eigenvalues  and  eigenvectors 
of  the  matrix  K-k**2D,  D=D(n)I 

c)  a(k,n)  are  constants  determined  by  initial  conditions. 

The  eigenvalues  and  eigenvectors  are  found  by  solving 

det ( K-k**2D- 1 1 ) =0  (8) 

The  ability  of  the  system  to  exhibit  pattern  formation 
depends  on  these  eigenvalues,  which  depend  in  turn  on  the 
rate  coefficients  in  K,  the  diffusion  coefficients  in  D,  and 
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the  "pattern  size  factor"  (or  spatial  eigenvalue)  K.  The 
eigenvalues  are  in  general  complex;  the  real  part  is  a 
growth  factor  and  the  imaginary  and  oscillatory  frequency 
associated  with  the  pattern.  If  any  eigenvalues  have 
positive  real  parts,  the  system  is  said  to  be  unstable  with 
respect  to  pattern  of  size  s=2/k,  because  any  disturbance  of 
the  steady  state  containing  an  element  of  this  will  grow 
exponentially.  If  all  eigenvalues  have  negative  real  part, 
the  system  is  stable. 

While  the  departure  from  equilibrium  remains  small 
enough  to  ensure  correctness  of  the  linearized  stability 
study,  the  pattern  size  for  which  the  growth  factor  is 
largest  will  dominate  the  disturbance.  Thereafter  nonlinear 
effects  come  into  play  (Gmitro  and  Scriven,  1966).  If  the 
eigenvalue  having  largest  real  part  is  entirely  real,  a 
"standing  wave"  of  concentration,  in  which  concentration  at 
each  point  steadily  increases,  will  be  seen;  if  it  is 
complex,  a  "travelling  wave"  will  be  seen,  in  which 
concentrations  at  each  point  cycle  with  frequency  given  by 
the  imaginary  part  and  grow  in  amplitude. 

In  summary,  a  stability  analysis  of  a  reaction- 
diffusion  system  shows  the  following  results.  A  spatially 
uniform  steady  state  of  the  system  is  subject  to  small 
disturbances,  whether  from  external  influence  or  internal 
"thermal  noise"  or  "molecular  fluctuations".  The  stability 
of  the  steady  state  in  response  to  such  disturbance  depends 
on  the  reaction  and  diffusion  coefficients  of  the  system;  if 
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these  are  such  that  instability  ensues,  spatial  patterns 
will  manifest  themselves.  The  nature  of  these  patterns 
depends  on  the  spatial  domain  of  the  system  and  the 
eigenvalues  of  K-k**2D. 

Gmitro  and  Scriven  (1966)  give  a  detailed  analysis  of 
the  patterns  which  may  arise  in  a  line,  a  circle,  a 
cylinder,  a  ring,  a  plane  (both  in  rectangular  and  polar 
coordinates)  and  a  sphere.  That  is,  they  give  the 
eigenvalues  or  elementary  patterns  for  the  Laplacian  for 
each  domain,  and  the  values  of  k**2  when,  in  a  closed  domain 
such  as  a  circle  with  conditions  imposed  on  the  boundary,  it 
must  assume  discrete  values  in  order  for  patterns  to  emerge. 
They  undertake  a  linearized  stability  analysis  for  one,  two 
or  three  participating  species  (corresponding  to  the  number 
of  coupled  reaction-diffusion  equations). 

3.3.6  Applications  of  Reaction-Diffusion  Kinetics 

A  brief  description  of  instances  in  which  reaction- 
diffusion  kinetics  have  been  used  in  modelling  pattern 
formation  is  given  below.  Some  of  these  models  will  be 
discussed  in  more  detail  when  computer -or i ented  models  are 
descr i bed . 

Turing  (1952)  investigates  the  departure  from 
equilibrium  of  linear  reaction-diffusion  systems  by  means  of 
linearized  stability  analysis.  Linear  systems  operating  in 
a  ring  of  cells  or  sphere  are  demonstrated  to  develop 
standing  or  travelling  waves  of  concentration  of 


. 

. 

' 


72 


exponentially  growing  amplitude;  a  discussion  of 

gastrulation  in  a  mammalian  embryo  is  linked  to  wave 
formation  in  a  sphere.  However  in  linear  systems  no 

limiting  of  the  exponential  growth  factor  is  possible  and 
thus  no  new  stable  state  or  pattern  can  be  established. 
Nonlinear  terms  coming  into  play  as  the  system  gets  further 
from  the  equilibrium  state  appear  to  balance  this 
exponential  growth  of  amplitude  in  the  wave  pattern  and 

allow  a  stable  pattern  to  be  established  (Gmitro  and 
Scriven,  1966;  Nicolis  and  Prigogine,  1977;  Grierer  and 

Mei nhar  t ,  1972). 

Crick  (1970)  has  investigated  the  feasibility  of 
reaction-diffusion  as  a  means  of  pattern  formation  in 
embryogenes i s .  Assuming  a  typical  size  of  a  primary  field 
as  0(10)  (between  about  30  and  100)  cells,  and  the 
establishment  of  the  field  (pattern)  within  5-10  hours,  he 
estimates  the  size  of  a  cell  and  rate  of  diffusion  necessary 
if  diffusion  is  to  be  assumed  the  means  of  spatial 

communication  of  intercellular  differences  in  concentration. 
He  obtains  a  value  for  the  diffusion  rate  which  is  of  the 
same  order  as  that  of  a  typical  metabolite  (c-AMP),  and  a 
cell  size  of  10-30  micrometres,  which  is  realistic.  The 
conclusion  is  therefore  that  the  parameters  of  the  diffusion 
process  are  such  that  it  is  a  practical  mechanism  for 
establishment  of  primary  fields. 

Grierer  and  Meinhart  (1972)  and  Meinhart  (1977)  use  two 
coupled,  nonlinear  reaction-diffusion  equations  for 
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activator  and  inhibitor  species  in  modelling  pattern 
formation  in  organisms  which  may  be  regarded  as  one¬ 
dimensional.  The  models  assume  initial  gradients  of  the 
density  of  cells  which  are  sources  of  the  two  chemical 
species.  Then,  if  a  react i on-di f f us i on  system  is  set  up  in 
which  the  activator  (which  catalyzes  both  itself  and  the 
inhibitor)  has  a  lower  rate  of  diffusion  than  the  inhibitor 
(which  inhibits  only  the  activator),  patterns  of 
concentration  of  the  two  species  are  established  with  an 
accentuated  peak  at  one  end  of  the  one-dimensional  array. 
It  is  suggested  that  such  patterns  could  act  as  prepatterns 
for  long-term  differentiation  of  additional  cells  into 
sources  of  the  species.  Computer  implementation  of  the 
model  (by  unspecified  numerical  techniques)  is  compared  to 
quantitative  data  from  hydra,  and  suggests  that  the  model 
can  account  for  certain  duplication  phenomena. 

Use  of  an  initial  gradient  of  density  of  source  cells 
(which  are  cells  differentiated  from  the  general  population) 
in  order  to  obtain  a  pre-pattern  for  further  cell 
differentiation  suggests  that  the  model  is  one  of  pattern 
" rei nforcement "  rather  than  pattern  "formation". 
Establishment  of  the  initial  density  gradient  has  not  been 
explained.  Less  initial  disparity  is  necessary  in  the  work 
of  Meinhart  (1977),  where  the  field  is  polarized  by  a  single 
source  of  activator  and  inhibitor  at  one  end.  Such  a  model 
is  used  to  discuss  segment  formation  in  insect  embryos. 

A  sophisticated  application  of  linearized  stability 


* 

. 

’ 


74 


analysis  of  a  two-species,  nonlinear  reaction-diffusion 
system  in  explaining  the  sequential  formation  and  shapes  of 
compartments  has  been  proposed  by  Kauffman  et  al.  (1978). 
In  their  model,  the  temporal  eigenvalues  of  the  linearized 
R-D  system  are  related  to  the  size  of  the  growing  disc, 
which  they  propose  to  be  approximately  elliptical.  The 
steady-state,  spatially  uniform  solution  of  the  system 
pertains  except  when  the  disc  attains  a  sequence  of  critical 
sizes.  At  these  sizes,  a  particular  eigenfunction  of  the 
Laplacian  (which  in  an  ellipse  are  the  Mathieu  functions) 
fits  exactly  into  the  ellipse  obeying  a  zero-flux  boundary 
condition;  thus  it  is  postulated  that  at  such  critical  sizes 
the  steady-state,  spatially  uniform  solution  becomes 
unstable  and  a  pattern  cor respondi ng  to  the  favoured 
eigenfunction  is  established.  Such  a  pattern  could  act  as  a 
prepattern  for  a  compar tmenta 1 i zat ion  event,  with  a 
compartment  boundary  specified  by  a  threshold  value  of  one 
of  the  species. 

As  the  disc  grows  larger,  the  pattern  no  longer  "fits" 
and  decays  to  the  uniform  situation,  until  the  next  critical 
size  is  reached  and  the  next  pattern  established.  The 
Mathieu  functions  (axes  of  the  ellipse,  confocal  hyperbolae, 
and  confocal  ellipses)  can  be  made  to  agree  in  general  shape 
and  timing  with  compar tmenta 1 i zat i on  in  the  wing  disc. 

Reaction-diffusion  relationships  of  a  primitive  Kind 
have  been  used  in  signalling  between  cells  in  computer  - 
oriented  models,  such  as  the  one  dimensional  model  of  Hydra 


V 


75 


developed  using  CELIA  (Herman  and  Rozenberg,  1975)  and  other 
models  based  on  L -systems  ( L i ndenmayer ,  1976).  Such 
modelling  of  intercellular  signalling  has  been  in  the  form 
of  uncoupled,  linear  difference  equations  and  will  be 
discussed  at  greater  length  in  the  next  section.  In 
general,  the  methods  fail  to  investigate  the  numerical 
stability  of  the  difference  equations  employed.  In 
consequence,  even  if  the  equations  were  of  sufficient 
complexity  to  admit  pattern  establishment  of  the  Kind  under 
discussion,  the  numerical  results  obtained  could  not 
demonstrate  them. 

Wi Iby  and  Ede  (1974)  have  developed  a  model  and 
simulation  for  differentiation  of  muscle  tissue  into  bone  in 
the  growing  chick  wing.  A  cell-to-cell  signalling  system  in 
terms  of  difference  equations  involving  one  species 
undergoing  synthesis,  "diffusion",  and  destruction  that  is 
dependent  on  two  thresholds  (for  initiation  and  termination) 
is  set  up  to  obtain  a  periodic  pattern  of  concentration  down 
the  limb.  Such  a  pattern  could  act  as  a  pre-pattern  for  the 
skeletal  pattern  of  the  limb.  The  authors  suggest  that  such 
a  periodic  distribution  is  more  plausible  than  a  situation 
in  which  a  monotonic  gradient  model  necessitates 
interpretation  of  a  series  of  thresholds  of  different 
magnitudes  being  i nterpretated  to  result  in  simular  cellular 
events  (formation  of  bone  cells).  The  difference  equations 
are  designed  to  produce  such  periodicity;  it  may  be  pointed 
out  that  the  eigenfunctions  of  the  Laplacian  on  a  rectangle 
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(corresponding  roughly  to  the  shape  of  the  limb  at 
intermediate  stages  of  development)  are  trigonometr ic  in 
form,  and  thus  use  of  a  reaction-diffusion  system  should 
allow  the  natural  establishment  of  a  periodic  pre-pattern. 

If  complicated  structures  such  as  the  limb  skeleton  are 
to  result  from  simple,  non-periodic  morphogen  gradients, 
then  one  mechanism  might  be  the  interaction  of  several  such 
gradients.  An  example  of  such  a  model  is  given  by 
MacWilliams  and  Papageorgiou  (1978),  who  introduced  a  pre¬ 
pattern  for  the  chick  wing  specified  by  thirteen  morphogens 
each  of  exponential  distribution.  Such  a  model  can  be 
constructed  from  the  fully  developed  limb  by  regarding  it  in 
the  light  of  a  two-  or  three-dimensional  map  of  contours 
around  the  various  structures.  The  shape  of  contours  around 
a  feature  determines  how  many  gradients  would  be  needed  to 
generate  such  a  field.  The  authors  suggest  that  cellular 
differentiation  at  such  a  feature  is  determined  by  a  "winner 
take  all"  rule  on  binding  the  morphogen  of  largest 
concentration  to  the  genome.  The  model  is  complex  and 
teleological  in  construction;  furthermore,  the  establishment 
of  the  numerous  gradients  is  not  discussed.  In  contrast,  a 
reaction-diffusion  system  with  few  interacting  species  is 
capable  of  establishing  equally  complex  patterns  in  a 
comprehens i b 1 e  manner . 

Grossberg  ( 1978a, b)  draws  an  analogy  between  the 
stability  character i st i cs  of  reaction-diffusion  systems 
( par t i cu 1 ar 1 y  those  involving  coupled,  long  term  inhibition 
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and  short  term  excitation  such  as  in  Grierer  and  Meinhart, 
(1972),  and  Kauffman  et  al.  (1978))  and  the  stability 
character i st i cs  of  shunting  networks  where  local  interaction 
in  governed  by  ordinary  differential  equations.  In  both 
types  of  systems  dynamically  changing  patterns  can  emerge 
from  initially  uniform  conditions.  Successful  approximation 
of  reaction-diffusion  systems  in  terms  of  finite  differences 
on  a  grid  produces  a  representat i on  mathematically 
equivalent  to  serial,  discrete  operation  of  mass-reaction 
laws  in  a  network  of  the  type  presented  by  Grossberg. 
Grossberg' s  approach  has  been  directed  towards  models  in 
population  biology,  ecology  and  psychophysiology,  with  some 
discussion  of  parallel  phenomena  in  developmental  biology. 
He  suggests  a  fundamental  class  of  control  mechanisms  in 
seemingly  competitive  systems  in  which  pattern  or  order 
appears  apparently  spontaneously. 

Various  investigations  have  been  made,  therefore,  of 
ways  in  which  biological  systems,  and  in  particular  primary 
fields,  can  develop  pattern  by  the  operation  of  small  scale, 
chemically  feasible  mechanisms.  The  success  with  which  such 
mechanisms  can  be  modelled  with  the  help  of  the  computer 
will  be  discussed  in  succeeding  chapters. 

3.4  Computer  Oriented  Models 

The  "real  system"  has  been  presented  as  a  mixture  of 
biology  and  existing  models  of  biological  phenomena.  A 
number  of  models  of  biological  system  in  which  simulation 
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plays  an  essential  part  have  been  constructed,  and  are 
described  below. 

3.4.1  Models  Based  on  L-Systems 

Developmental  algorithms  for  multicellular  organisms 
have  been  developed  based  on  the  work  of  Lindenmayer  on  L- 
systems  (Lindenmayer,  1975).  An  L-system  is  a  finite 
automaton  in  which  components  of  the  array  can  be  replaced 
by  arrays.  Because  of  the  complexity  of  the  geometry  of 
such  replacements  in  two  and  three  dimensions,  most  of  the 
successes  of  L-systems  have  been  in  the  modelling  of 
structures  of  a  one  dimensional  character.  Extension  to 
higher  dimensions  by  means  of  "graph-L-systems" ,  where 
components  are  nodes  which  can  be  replaced  by  nodes  and 
links,  is  being  attempted  (Lindenmayer,  1975). 

An  L-system  is  a  di screte/di screte  model:  components 
(cells)  are  discrete,  and  the  state  set  is  discrete.  The 
state  of  a  cell  usually  refers  to  a  genetic  state  which  is 
heritable  and  may  change  according  to  a  transition  function 
which  is  uniform  over  the  array,  the  present  state  of  the 
cell,  and  the  states  of  its  neighbours.  The  transition 
function  may  be  deterministic  or  non-determi ni s t i c . 

One-dimensional  L-systems  correspond  exactly  to  formal 
language  theory:  the  set  of  states  is  an  alphabet,  and  all 
developmental  possibilities  may  be  generated  by  the 
operation  of  the  transition  function  (set  of  productions)  on 
the  starting  state  (axiom).  Because  of  this  formal 
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equivalence,  many  predictions  about  L-systems  can  be  made 
theoretically.  Systems  are  classified  according  to  the 
connectivity,  the  number  of  neighbours  affecting  each  cell, 
and  have  been  used  to  model  appropriate  developmental 
systems  of  a  filamental  nature. 

The  work  of  Herman  et  al.  (1974)  applies  the 
techniques  of  L-systems  to  development  in  simple  two- 
dimensional  organisms.  Their  FORTRAN  program  CELIA 
simulates  linear  iterative  networks  of  identical  modules  (or 
cells).  Cells  have  a  finite  number  of  possible  states  and 
change  state  depending  on  their  own  previous  state  and 
threshold  values  of  signals  from  their  two  neighbours.  A 
complete  set  of  rules  for  change  of  state  must  be  supplied 
by  the  user.  The  end  cells  in  a  line  can  communicate  with 
the  environment.  Cells  divide  synchronously.  Branches  can 
be  grown  on  linear  organisms  by  means  of  a  bracketted  list 
s t ructure . 

CELIA  has  proved  useful  in  demonstrating  the 
effectiveness  of  a  linear  gradient  in  producing  coordinated 
changes  of  state  in  simulations  of  certain  organisms.  The 
developmental  rules  for  red  algae,  patterns  on  sea  snail 
shells,  and  grafting  experiments  in  Hydra  are  some  of  the 
areas  to  which  the  program  has  been  applied.  It  has  been 
demonstrated  that  a  polar  regulation  of  cell  states  can  be 
caused  without  imposing  polarity  in  individual  cells. 

The  limitations  of  CELIA  are  apparent.  Only  linear 
gradients  parallel  to  the  linear  iterative  arrays  can  be 
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investigated.  The  authors  themselves  point  out  the  need  for 
probabilistic  changes  of  state  and  cell  division.  Their 
simulation  is  implemented  in  FORTRAN,  making  their  data 
structures  difficult  to  maintain.  A  departure  from  strict 
synchronous  processing  in  such  structures  is  introduced  by 
Hogeweg  (1978),  in  which  advantage  is  taken  of  the  low 
activity  level  in  the  array  by  using  a  discrete  event 
technique  and  SIMULA  in  programming  the  model. 

3.4.2  Growth  Models  in  Arrays  of  Cells 

Wilby  and  Ede  (1976)  in  further  study  of  the  chick 
wing,  have  constructed  a  simulation  of  growth.  A 
rectangular  array  with  initial  cells  on  one  border  is 
subjected  to  various  gradients  of  mitotic  capability.  A 
cell  divides  probabilistically  if  its  internal  situation 
places  it  above  a  threshold  of  mitotic  ability  determined  by 
its  position.  The  shortest  path  to  an  unoccupied  position 
is  then  found,  and  all  cells  along  that  path  displaced 
outwards.  By  these  means  the  shapes  of  normal  and  abnormal 
wings  have  been  obtained,  showing  that  limb  shape  is 
possibly  mediated  only  by  controlled  cell  division. 

The  Wilby  and  Ede  simulation  is  an  elegantly  simple 
model  of  growth,  but  i nterce 1 1 u 1 ar  communication  and  the 
data  structures  needed  to  describe  it  are  not  considered. 

A  simulation  of  embryonic  development  by  Raven  and 
Bezem  (1976)  attempts  to  characterize  differences  between 
cells.  They  represent  cells  as  spheres  of  different  size, 


' 


. 

■ 


81 


with  cell  flattening  along  contact  surfaces  represented  by 
the  chord  of  intersection  of  two  spheres.  This  method 
permits  a  close  view  of  the  early  stages  of  development, 
when  there  are  not  many  cells.  By  assuming  that  the  type  of 
a  cell  depends  on  concentration  gradients  of  axial  and 
radial  shape,  they  produce  realistic  cell  differentiation 
patterns  in  a  simulation  of  Lymnaea. 

A  di screte/di screte  model  of  cell  division  which  is 
used  to  study  clonal  patterns  has  been  developed  by  Ransom 
(Ransom,  1977)  and  programmed  in  an  ALGOL-based  language 
Edinburgh  IMP.  Cells  do  not  change  state  (therefore  the 
model  is  not  truly  developmental)  but  are  defined  by  a 
heritable  clonal  marker  and  a  position  in  an  array  of  fixed 
maximum  size.  The  array  is  hexagonal ly  arranged.  Cells  are 
processed  sequentially  in  each  time  step;  a  cell  is  allowed 
to  divide  if  it  has  not  divided  already  in  the  current  time 
step,  and  if  its  division  does  not  cause  separation  of  any 
pair  of  cells  having  the  same  clonal  marker  ("stickiness" 
rule).  When  a  cell  divides,  the  nearest  free  space  is  found 
and  cells  along  that  direction  are  displaced  outward.  The 
displacement  is  simulated  by  moving  the  entries  for  cells  in 
a  hash  table  representing  the  array. 

Cells  near  the  centre  of  a  population  are  less  likely 
to  divide  in  this  model  as  they  are  more  likely  to  violate 
the  "stickiness  rule"  than  cells  near  the  edge.  Ransom  has 
used  his  model  to  simulate  growth  in  Drosophila  imaginal 
discs,  in  which  boundary  cells  do  divide  more  frequently 


. 

’ 


82 


(Ransom,  1977b).  Division  in  the  direction  of  nearest  free 
space  implies  division  in  radial  directions  in  a  "circular" 
disc,  and  Ransom  argues  that  straight  radially  oriented 
clonal  boundaries  observable  at  some  stages  of  disc 
development  may  be  produced  by  such  directed  cell  divisions. 

Ransom's  model  is  not  genetically  oriented;  cells  do 
not  change  state,  nor  is  their  state  described  by  more  than 
one  discriminant.  Furthermore  they  do  not  communicate. 
Moreover  cells  divide  without  respect  to  local  conditions  or 
their  own  state.  The  only  factor  affecting  division  is  the 
"stickiness  rule",  in  which  an  event  remote  from  a  cell  may 
prevent  its  division. 

Duchting  (Duchting,  1977)  has  constructed  a  model  of 
growth  in  a  10x10  matrix  of  cells  using  the  mi croprocessor 
system  INTEL  8080.  It  bears  some  similarity  to  Conway's 
Game  of  Life  in  that  cells  become  active,  die  or  divide 
according  to  the  activity  status  of  their  neighbours.  Cells 
inherit  a  marker  which  also  indicates  a  division  rate.  The 
model  has  been  used  to  simulate  cell  populations  of 
different  division  rates  in  competition,  and  may  therefore 
be  of  utility  in  studying  the  prol i ferat ion  of  cancer  cells. 

3.4.3  Intracellular  Growth  Models 

Various  models  have  been  constructed  which  attempt  to 
describe  the  internal  factors  of  a  cell's  biochemistry  which 
determine  its  division  cycle  time.  Simulation  of  such 
models  results  in  the  collection  of  statistics  about  the 
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cells  in  a  population  such  as  the  distribution  of  division 
times,  and  are  not  spatially  oriented.  Typical  of  such 
models  is  that  of  Alberghina  (Alberghina,  1977),  in  which  a 
cell  divides  if  its  production  of  ribosomes  and  proteins 
reaches  a  threshold  level.  The  environment  can  affect  the 
cell  by  adjustment  of  reaction  rates  in  the  equations 
modelling  the  biochemical  processes. 

While  i ntrace 1 1 u 1 ar  models  of  cell  division  could  be 
used  in  looking  at  a  population  of  cells  as  a  whole,  it  is 
more  likely  that  statistics  generated  by  such  models  would 
be  used  directly  in  obtaining  a  probaba  1 i s t i c  model  of 
growth  in  a  population,  thus  "lumping"  the  biochemical 
mechanisms  within  the  cell.  Such  techniques  are  employed  in 
Chapter  Six. 

3.5  The  Experimental  Frame 

The  task  of  recognizing  the  experimental  frame,  through 
which  the  real  system  is  to  be  viewed  for  modelling  and 
simulation,  is  that  of  identifying  both  the  information 
going  into  populations  of  cells  in  a  developing  organism  and 
the  measurable  results  which  exper i menter s  find  in  such 
organisms.  Such  input  and  output  will  be  peculiar  to  the 
particular  system  under  study,  but  some  general  remarks  may 
be  made  concerning  experimental  frames  for  developmental 
systems  of  interest. 

In  order  that  organised  patterns  of  differentiating 
cells  may  be  visualized  with  clarity,  an  initial  decision  is 
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made  to  study  systems  which  are  either  inherently  two- 
dimensional  or  may  be  viewed  in  two-dimensional  cross- 
section.  A  drosophila  disc  is  essentially  a  two-dimensional 
structure,  while  a  chick  wing's  pattern  of  cartilage  and 
muscle  may  be  regarded  in  sections  perpendi cu 1 ar  to  the 
dorsa 1 /vent r a  1  axis  or  to  the  proximo/distal  axis.  It 
should  be  noted  that  a  base  model  of  development  should  not 
reflect  this  limitation  to  two  dimensions,  though  lumped 
models  constructed  within  experimental  frames  will  do  so. 

The  most  significant  limitation  of  the  experimental 
frame  arises  from  developmental  geneticists'  use  of  "top- 
down"  or  difference  methods  in  studying  developing  systems. 

A  perturbation  or  sequence  of  naturally  occurring  events 
stimulates  division  of  an  initial  configuration  into 
subpopulations  of  cells  with  mutual  differences.  Those 
cellular  activities  which  are  common  to  all  cells,  and 
remain  common  to  all  cells  during  the  experiment,  are  not 
character i sed .  The  experimental  frame  implicit  in  such 
approaches  neglects  all  inputs  and  char acter i s t i cs  of  cells 
which  are  not  relevant  to  the  developmental  changes  under 
study. 

In  general,  inputs  to  systems  consist,  then,  of 
"perturbations"  (for  example,  an  external  chemical  signal  to 
a  population  of  cells,  or  a  surgical  cut  or  other 
disturbance  within  a  population).  Inputs  to  individual 
cells  will  be  signals  thought  to  be  involved  in  the 
formation  of  patterns.  Outputs  may  be  of  a  qualitative 
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nature,  such  as  the  shapes  of  developing  systems  or  the 
normality  of  an  emerging  pattern.  Further  discussion  of 
experimental  frames  relevant  to  particular  systems  will  be 
undertaken  in  Chapter  Six. 
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CHAPTER  FOUR 


The  Modelling  Framework 


While  a  particular  developmental  system  of  the  type 
discussed  in  the  previous  chapter  will  be  char acter i zed  by  a 
particular  model,  there  are  aspects  of  such  models  which  are 
held  in  common.  In  particular  the  base  model  of  the  genetic 
mechanisms  operating  in  cells  is  almost  certainly  universal, 
to  the  extent  that  inherited  genetic  information  and 
environmental  factors  interact.  A  particular  model  for  a 
particular  system  exists  within  a  modelling  framework, 
sharing  some  elements  with  other  systems  and  having  some 
unique  characteristics.  The  modelling  framework  as  a  whole 
may  be  looked  upon  as  a  pool  of  models,  or  collection  of 
options,  correspondi ng  to  various  aspects  of  the 


developmental  process,  different  developmental  systems, 
alternative  explanations  of  developmental  processes. 
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4.1  The  Underlying  Base  Model 

Consider  a  cell  in  a  population  of  cells.  It  has 
inherited  state  of  determination  for  the  activation  of 
portions  of  its  DNA.  This  state  of  the  DNA  is  indicative  of 
the  cell's  present  and  future  chemical  activity.  The 
present  expression  of  the  DNA  will  be  referred  to  as  the 
cell's  'functional  nature'.  Two  cells  which  are  in 
different  determi nat i ona 1  states  (have  different  activation 
of  DNA)  have  at  least  potentially  different  functional 
natures . 

Consider  a  population  of  cells  which  are  mutually 
undifferentiated.  They  are  in  the  same  state  of  activation 
of  their  DNA,  and  they  have  the  same  functional  nature.  Now 
an  event  occurs  which  divides  the  population  into  two  groups 
with  different  activation  states  in  their  DNA;  that  is,  they 
are  distinct  in  their  determination.  It  may  not  be  Known 
what  the  difference  between  them  is,  but  at  least  one  such 
difference  exists.  The  difference,  whether  it  be  one  or  a 
hundred  sites  of  different  activation  on  the  DNA,  is  the 
result  of  one  determi nat i ona 1  event  and  can  be  symbolized  by 
one  "switch”.  Such  a  switch  may  be  likened  to  a  sensor  gene 
in  the  Britten  and  Davidson  model.  The  sensor  gene  acts  as 
a  Key  to  the  activation  of  other  levels  of  control  over  a 
certain  number  of  cell  products. 

The  situation  is  pictured  in  Fig.  4.1.  Let  the  two 
groups  of  cells  be  A  and  B,  and  assign  switch  codes  of  1  0 
and  0  1  to  them  respectively.  The  fact  that  the  codes  are 
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FIGURE  4.1  Coding  in  Subdividing  Populations 
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different  reflects  the  occurrence  of  a  determi nat i ona 1 
event.  The  1  in  each  code  reflects  the  activation,  or 
enabling  for  future  activation,  of  specialized  functions  for 
that  group  of  cells.  (If  the  subdivision  had  been  tertiary, 
codes  of  1  0  0,  0  1  0,  and  0  0  1  could  have  been  assigned.) 
The  Hamming  distance  between  the  codes  for  A  and  B  cells  is 
2;  this  reflects  the  fact  that  2  things  would  have  to  take 
place  if  A  cells  were  to  become  B  cells:  a  return  to 
indeterminate  state  ( 1  0  to  0  0 )  and  re-determi nat i on ( 0  0  to 
01).  In  terms  of  a  representat i on  of  determination  as  a 
decision  tree,  the  "distance"  between  two  groups  of  cells  at 
a  particular  moment  in  time  (at  the  same  level  in  the  tree) 
may  be  regarded  as  the  number  of  branches  in  the  shortest 
path  connecting  them. 

It  may  be  noted  that  the  codes  1  0  and  0  1  may  model 
equally  well  a  repression  of  activity  as  indicated  by  the  0 
in  each  code.  In  a  repression  model,  the  indeterminate 
state  would  be  represented  by  1  1 . 

Similarities  between  this  character izat ion  and  that  of 
Kauffman  (1973)  will  be  noted;  however,  the  code  proposed 
here  is  intended  to  have  greater  significance.  Kauffmann 
assigned  his  codes  at  random,  the  only  condition  being  that 
Hamming  distance  should  correlate  with  frequency  data.  In 
the  proposed  model  the  code  for  a  group  of  cells  accumulates 
as  time  goes  on  and  determi nat i ona 1  events  occur.  The  code 
is  thus  a  history  of  each  cell  in  reference  to  all  other 
cells  since  they  were  relatively  indeterminate. 
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Furthermore,  the  code  for  a  group  of  cells  is  a  mapping  onto 
the  DNA ,  as  each  1  represents  activated  or  enabled  portions 
of  the  DNA. 

Assume  for  example  that  the  different  functional  nature 
of  the  cell  groups  is  to  at  least  some  degree  expressed  in 
terms  of  present  functional  activity.  It  may  not  be  Known 
what  precise  chemical  activity  is  exclusive  to  cells  of  type 
A  (1  0),  but  what  is  Known  is  that  all  A  cells  have  the  same 
functional  nature  prior  to  any  subsequent  determi nat i ona 1 
event.  That  is,  all  A  cells  are  doing  (or  can  do)  the  same 
things,  while  cells  of  type  B  are  doing  (or  can  do)  at  least 
one  thing  that  is  different.  It  is  easy  to  see  that 
gradients  of  information  can  be  set  up  as  a  result  of  this 
differential  activity  of  cells  in  the  population.  Consider, 
for  example,  diffusion  through  the  entire  population  of  a 
chemical  produced  only  by  cells  of  type  A.  A  simulation 
experiment  could  be  set  up  to  see  whether  the  pattern  of 
concentration  resulting  from  such  a  simple  diffusion 
exhibited  an  anomaly  along  the  line  W V  in  Fig.  4.1.  Thus, 
without  Knowing  the  actual  nature  of  activity  of  cells  of 
type  A,  it  is  possible  to  study  patterns  arising  from  those 
cells'  expression  of  their  nature  in  the  form  of  signals  to 
the  entire  population  or  portions  thereof. 

If  clues  do  exist  as  to  the  functional  nature  of  cells 
of  type  A,  this  information  can  be  used  in  determining  how 
signals  would  be  transmitted.  For  example,  simple  diffusion 
might  be  appropriate.  If  data  are  available  about  the 
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nature  of  boundaries  between  compartments,  such  information 
will  again  be  reflected  in  the  way  information  is  allowed  to 
travel.  If  it  is  Known  what  shape  divisions  will  arise,  and 
when,  different  transmission  possibilities  may  be  tried  in 
order  to  throw  light  on  what  the  actual  signals  may  be.  In 
this  way,  a  model  and  experimental  data  can  be  mutually 
sens i t i ve . 

Each  cell  in  a  population  has  the  above  code  of  binary 
digits  which  symbolizes  its  state  of  determination  with 
respect  to  the  rest  of  the  population.  It  shares  this  code 
with  all  cells  like  itself.  In  addition,  however,  it  has 
local  information  in  the  form  of  concentrations  of 
chemicals,  fields  of  ionic  potential,  pressure,  and  other 
physical  realities  in  the  cell's  environment.  These 
quantities  may  be  loosely  classified  into  3  groups: 

(1)  external  forces,  e.g.  gravity,  organism-wide  hormonal 
flow; 

(2)  processes  common  to  all  cells,  e.g.  fundamental  life 
support  metabolic  functions; 

(3)  local  levels  of  specialized  functions;  a  cell  may  have 
its  code  in  common  with  certain  other  cells  and  yet  differ 
from  them  in  the  concentration  of  products  activated  by  the 
code,  perhaps  because  of  positional  differences  or  different 
timings  of  the  onset  of  new  functional  activities. 

Using  the  proposed  difference  model,  as  many  of  the 
physical  influences  in  the  above  three  groups  as  are  thought 
to  be  of  interest  can  be  explicitly  considered  in  designing 
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a  particular  modelling  experiment;  that  is,  an  appropriate 
experimental  frame  can  be  selected.  For  each  of  these 
chosen  quantities  an  account  of  changes  occurring  can  be 
Kept  as  a  state  parameter  of  the  cell.  The  code  of  digits 
will  be  a  Key  to  which  of  the  functions  of  type  (3)  a 
particular  cell  is  pursuing  actively  (it  may  be  receiving 
signals  from  cells  with  different  codes).  The  concept  of 
each  cell's  state  can  be  illustrated  with  two  examples. 

(1)  A  population  of  cells  of  uniform  determination  (no 
code)  is  subjected  to  the  directed  flow  of  a  hormone,  and  a 
pressure  gradient  exists,  due  to  gravity,  which  influences 
this  flow. 

Each  cell  is  char acter i zed  by  its  concentration  of  the 
hormone  and  the  local  pressure.  The  value  of  these  as  a 
function  of  time  is  given  by  the  solution  of  an  equation  of 
hydrodynamic  flow  across  the  population.  While  other 
processes  are  in  operation  in  the  cells,  these  two  are  the 
focus  of  attention.  The  desire  may  be  to  investigate 
patterns  arising  in  the  concentration  of  hormone. 

(2)  A  disc  with  two  compartments  is  viewed  in  terms  of 
the  differential  functional  activity  present. 

Production  of  A-specific  signals,  Keyed  by  the  1  in  the 
code  1  0,  starts  at  the  determination  event,  while  B 
signals,  Keyed  by  the  1  in  the  code  0  1,  start 
simultaneously.  The  O' s  indicate  input-only  fields  for 
signals.  If  each  cell  produces  a  unit  quantity  of  signal  in 
a  unit  of  time,  patterns  of  relative  concentrations  of 
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signal  can  be  studied.  Each  cell  is  characterized  by  the 
concentration  of  the  two  signals. 

A  genetic  model  has  been  proposed  here  which  views  a 
developmental  system  as  a  spectrum  of  cells,  each  type  of 
which  displays  a  certain  behavior  reflecting  its 
determi nat i ona 1  state,  the  pattern  of  activation  in  its  DNA. 
This  model  is  the  base  model  for  developmental  systems.  The 
base  model  is  designed  to  embody  the  ways  in  which  the 
activational  state,  as  symbolized  by  the  code  of  genetic 
switches,  interacts  with  the  behavior  to  produce 
determi nat i ona 1  events  in  an  organized  sequence.  The  model 
is  dynamic,  flexible  in  the  scope  of  its  attention  (the 
whole  organism  or  any  subgroup  of  cells),  and  quantitative. 

A  model  proposed  by  Garci a- Be  1 1 i do  (1975)  bears  an 
interesting  similarity  to  the  ideas  presented  above.  As  a 
result  of  his  experimental  work,  he  has  identified  the 
"homeotic"  genes  as  controlling  the  development  of  whole 
compartments;  that  is,  he  has  proved  the  existence  of  a 
single  gene  controlling  the  action  of  a  battery  of  producer 
or  "realizator"  genes  in  a  manner  similar  to  that  envisaged 
by  Britten  and  Davidson  (1969).  The  homeotic  gene,  he 
suggests,  is  switched  on  by  the  product  of  an  "activator" 
gene  in  conjunction  with  an  "inducer"  molecule.  The 
activator  gene  is  presented  as  a  means  of  reflecting  the 
genetic  inheritance  or  determined  state  of  the  cells,  while 
the  inducer  molecule  is  a  means  of  involving  local 
information,  thus  obtaining  positional  specificity  of  the 
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homeotic  gene.  The  inducer  could  be  a  signal  from  outside 
the  population  of  cells,  or  a  signal  from  cells  of  a 
different  type.  The  quantitative  model  suggested  above  uses 
the  present  activation  code  to  show  the  determined  state  of 
the  cells,  while  local  information,  as  achieved  through  the 
functional  quantities  maintained  for  each  cell,  is  used  to 
obtain  positional  specificity  of  new  activation  code.  In 
essence,  the  two  models  have  the  same  goal,  to  show  the 
interaction  of  present  DNA  activation  and  local  information 
in  producing  new  DNA  activation. 

In  further  support  of  the  approach  taken,  it  may  be 
noted  that  Lewis  (1976)  has  discussed  the  discrete  nature  of 
compar tmenta 1 i zat i on  events  in  response  to  threshold  values 
of  continuous  gradients  of  information.  He  concludes  that 
such  development  is  appropr i ate  1 y  modelled  in  terms  of 
"genetic  switching  circuits". 

Restriction  of  the  base  model  by  the  experimental  frame 
results  in  char acter i zat i on  of  cells  and  cell  populations  in 
terms  of  a  difference  model;  that  is,  only  those  portions  of 
the  genetic  code,  and  those  related  activities,  which 
express  accumulating  differences  between  cells  during  a 
developmental  process,  are  maintained  explicitly. 

Developmental  events  occur  in  an  organism  because  of 
the  recognition  of  the  attainment  of  critical  values  of 
certain  behavioral  parameters  in  the  cell  population. 
Recognition  of  such  an  event  occurs  when  critical  values  are 
reached  in  a  particular  subcollection  of  cells,  along  a 
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curve  which  becomes  a  compartment  boundary,  for  example. 


4.2  Common  Elements  of  Lumped  Models 

There  are  two  significant  simplifications  common  to  all 
cases  of  the  production  of  lumped  models  capable  of 
translation  into  iterative  specifications  i nterpretab 1 e  by 
digital  computers.  First,  any  continuous  time  base  has  to 
be  approximated  by  a  discrete  one.  The  appropriate  length 
of  the  time  step  will  depend  on  the  particular  model  being 
constructed.  Second,  in  the  absence  of  true  parallel 
processing  for  individual  cells,  any  activities  which  are  in 


reality  simultaneous 

mus  t 

be 

approx i mated 

by  serial 

processing.  The 

degree 

to 

which 

this  approximation  is 

successf u 1  will 

depend 

on 

the 

treatment  of 

par  t i cu 1 ar 

processes  in  particular  models.  The  validity  of  these 
simplification  steps  must  be  decided  in  behavioral  terms; 
that  is,  in  terms  of  the  relationship  between  a  particular 
lumped  model,  in  a  particular  experimental  frame,  to  the 
base  model  and  real  system. 

The  fundamental  genetic  mechanisms  in  the  base  model  of 
development  translate  into  structures  which  are  common  to 
all  1 umped  mode  1 s . 

4.2.1  Characterization  of  Individual  Cells 

A  cell  is  an  entity  which  has  local  information  of  both 
quantitative  and  operational  Kinds.  At  a  given  time,  the 
state  of  the  cell  is  determined  by  certain  local  values  of 


96 


variables  including  global  forces  such  as  pressure, 
metabolic  levels,  and  concentrations  or  "signal  levels" 
cor respondi ng  to  expression  of  the  genetic  nature  of  the 
cell  and  its  neighbours.  In  a  particular  modelling 
situation,  a  subset  of  these  variables,  those  directly 
involved  in  the  developmental  process  of  interest,  will  be 
recognized  as  the  state  set  of  the  cell. 

A  cell  receives  input  from  its  immediate  environment  in 
the  form  of  information  about  the  state  of  those  neighbours 
with  which  it  is  in  contact.  If  it  is  in  physical  proximity 
to  another  cell,  but  a  barrier  to  certain  signal  exists 
between  them,  no  input  of  these  values  is  transmitted.  The 
cell  communicates  information  about  its  own  state  to  its 
neighbours  in  a  similar  manner. 

A  cell  contains  operational  information  in  the  form  of 
its  activation  code.  This  information  is  a  set  of  switches 
indicating  which  of  its  variables  change  because  they  form 
an  active  part  of  the  cell's  expression  of  its  genetic 
nature.  These  indicated  variables  are  manipulated  by  the 
cell  itself  and  do  not  merely  change  passively  due  to 
signalling  between  cells. 
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4.2.2  Characterization  of  a  Group  of  Cells  with  the  Same 
Genetic  Code 

A  group  of  cells  of  identical  genetic  nature  (having 
the  same  set  of  switches)  comprises  an  element  in  the 
formation  of  patterns  in  the  developing  organism.  The 
group's  component  cells,  therefore,  are  considered  linked 
together  in  a  structure  which  has  shape  and  position  in  the 
total  population.  A  "code  group",  then,  is  defined  by  a 
list  of  cells  of  identical  genetic  state,  together  with 
their  positions,  and  shape  parameters  of  the  whole  group. 
Description  of  shape  is  dependent  on  the  particular  system; 
for  example,  if  the  group  is  a  coherent  patch,  its  shape 
could  be  described  by  a  list  of  perimeter  cells  and  their 
pos i tons . 

A  code  group  may  itself  be  regarded  as  a  system  with  a 
set  of  state  variables,  as  a  single  cell  is  defined  at  a 
lower  level  in  the  modelling  hierarchy.  One  element  of  the 
state  variable  set  is  the  string  of  genetic  switches.  The 
genetic  activation  code  changes  upon  the  recognition  of  a 
developmental  event  in  the  component  cells  of  the  group. 
Another  element  of  the  state  set  is  the  shape  description. 

A  perimeter  list  will  change  during  growth  of  the  group  due 
to  cell  division,  or  because  of  division  of  the  group  into 
subgroups  during  an  event  such  as  compar tmenta 1 i zat i on . 

Thus  changes  in  the  state  variables  of  a  group  are 
mediated  by  events  on  the  cellular  level.  It  is  appropriate 
in  most  modelling  situations  to  define  a  code  group  as 
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autonomous  (input  free)  and  to  consider  input  only  at  the 
cellular  level.  Output  of  a  code  group  would  depend  on  the 
particular  model;  it  might  be,  for  example,  the  shape 
descr i pt ion . 

4.2.3  Char acter i za t i on  of  a  Developmental  System 

A  developmental  system  is  a  collection  of  code 
groups  of  cells,  each  group  having  a  distinct  genetic 
nature.  The  system  may  be  described  by  listing  these  groups 
and  locating  them  relative  to  each  other.  The  shape 
description  contained  within  a  group  may  then  be  used  to 
construct  a  picture  of  the  developmental  system  as  a  whole. 

A  developmental  system  changes  by  the  addition  of  new 
groups  (cor respond i ng  to  newly  specialized  cell  types)  to 
the  list  (or,  conceivably,  by  deletion  of  groups). 

4.2.4  Cell  Di vi sion 

Growth  is  a  process  inseparable  from  development.  The 
formation  of  patterns  of  cell  types  in  an  organism  depends 
on  the  mechanism  of  cell  pro! i ferat ion  as  well  as  on  events 
of  a  genetic  nature  (Wolpert,  1978).  Any  spatial 
configuration  must  depend  on  shape  and  the  rate  of  expansion 
of  the  organism.  Modelling  the  growth  process  is  therefore 
an  essential  part  of  modelling  any  developmental  system. 

The  shape  of  an  organism  is  controlled  by  the  location, 
timing  and  direction  of  cell  divisions.  A  particular  cell's 
decision  to  divide  is  made  on  the  basis  of  an  intracellular 
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model  (for  example,  a  cell  cycle  clock,  or  threshold  value 
of  a  metabolite)  or  on  the  basis  of  an  i nterce 1 1 u 1 ar  model 
(for  example,  threshold  level  of  a  "growth  hormone",  or  a 
probabilistic  model  of  division).  If  a  cell  does  decide  to 
divide,  it  may  do  so  in  a  preferred  direction,  for  example 
in  the  direction  of  least  pressure,  or  down  some  other 
gradient.  A  pool  of  growth  models  containing  such  options 
must  be  made  available  to  aid  in  the  construction  of  models 
of  particular  developmental  systems. 


4.3  Formal  Characterization  of  Lumped  Models 
4.3.1  Specification  of  Individual  Cells 
A  cell  is  a  structure 


C=<T,X,W,Q,f> 

where 

T  is  the  time  base 

X  is  the  input  value  set  (range  of  possible  values  of 
signals  entering  the  cell) 

W,  a  subset  of  (X,T),  is  the  input  segment  set 
Q  is  the  state  variable  set 

f  is  the  state  transition  function  f:QxW  -->  Q 

The  function  f  yields  the  state  at  the  end  of  each 
input  segment.  C  is  an  input-output  system  if  W  is  closed 
under  composition  (successive  I/O  experiments  can  be 
performed)  and  f  has  the  composition  property,  that  is,  for 
w(1),w(2)  6  W  and  contiguous, 

f(q,w( 1 ) »w ( 2 ) )=f(f(q,w(1 ) ) , w ( 2 ) ) 
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Examp  1 e : 

T  is  a  discrete  time  base 

X  =  { z ( t ) | 0  <  z(t)  <  Z } ,  the  concentration  of  substance  x 
entering  the  cell  from  its  neighbours  in  time  step  t 
W  is  a  subset  of  (X,T) 

Q  is  the  state  variable  set,  Q=x,  the  concentration  of 
substance  x  in  the  cell 

Y  is  the  output  value  set,  { y | y  an  integer,  x-1<y<x} 

f  is  the  state  transition  function,  modelling  the 
concentration  changes  due  to  reaction  and  diffusion  in  the 
system. 

4.3.2  Specification  of  Code  Groups 
A  code  group  is  a  structure 
G=<T , Q , Y , f  , g>  where 
T  is  a  time  base 
Q  is  the  state  variable  set 

Y  is  the  output  value  set 

f  is  the  state  transition  function,  f:Q  -->  Q 
g  is  the  output  function  g:Q  -->  Y 
Example:  Let 

T  be  a  discrete  time  base 

Q  =  (s,L,D)  where  s  is  a  string  of  binary  digits  symbolizing 
the  genetic  activation  state  of  the  group; 

L  is  the  list  of  cells  in  the  group  with  their  addresses; 

D  is  the  list  of  perimeter  cells. 

Suppose  two  Kinds  of  processes  are  going  on  in  the 


p 

. 

. 
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group,  cell  division  and  production  of  a  substance  x. 
Suppose  that  if  concentration  of  x  exceeds  a  critical  level 
c  in  over  80%  of  the  cells  a  genetic  event  will  occur 
changing  the  value  of  s  from  s ( 1 )  to  s(2).  Then  f  may  be 
given  by 

f(s(1),L(1),D(1))=(s(2),L(1),D(1))  if  the  genetic  event 

occurs 

=  ( s ( 1 ) , L ( 2 ) , D ( 2 )  )  otherwise 

where  L ( 2 )  is  the  list  of  cells  and  daughter  cells  and  D ( 2 ) 
is  the  new  perimeter  list  after  cell  divisions  during  time 
t . 

4.3.3  Specification  of  a  Developmental  System 
A  developmental  system  is  a  structure 
<T  ,  Q  ,  f >  where 
T  is  a  time  base 
Q  is  the  state  variable  set 

f  is  the  state  transition  function  f:Q  -->  Q 
Example:  Let 

T  be  a  discrete  time  base; 

0  =  L,  a  list  of  code  groups  with  associated  location 
parameters . 

Suppose  L  changes  whenever  a  code  group  undergoes 
subdivision  into  two  code  groups.  Then  f ( L ( 1 ) )  =  L ( 2 ) 
where  L ( 2 )  is  L ( 1 )  minus  the  dividing  group  plus  the 
two  new  subgroups . 
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f(L(1))=L(1)  if  no  genetic  event  occurs  in  t 
=L(2)  otherwise. 

4.4  Variable  Elements  For  Individual  Models 

It  may  be  seen  that  while  the  structure  of  the 
cell,  the  code  group  and  the  developmental  system  are 
determined  by  the  modelling  framework,  the  individual 
elements  within  the  structure  are  not  constrained.  In 
the  cell,'  all  the  sets  T,X,W,Q,Y  are  free  to  be  defined 
for  a  particular  developmental  model.  The  transition 
function  f  is  constrained  only  by  considerations  such 
as  a  desire  that  it  obey  the  composition  property  (if  a 
cell  is  to  be  an  I/O  system)  or  to  be  time  invariant. 
Similarly  the  elements  of  the  code  group  and  organism 
structures  may  be  defined  freely. 

The  formal  specification  of  these  structures 
cannot  be  given  in  iterative  form  because  of  the  lack 
of  specification  of  these  free  elements,  in  particular 
of  the  time  base.  When  experiments  are  performed  on 
particular  models,  a  formal  iterative  specification 
with  all  free  elements  fixed  will  be  given. 
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CHAPTER  FIVE 

Implementat ion 

Corresponding  to  the  modelling  framework  describing 
structures  and  elements  common  to  all  the  particular  models 
of  developing  systems  to  be  constructed,  there  is  a 
computing  framework  implementing  these  common  features.  The 
framework  consists  of  modules  incorporating  the  separate 
entities  and  processes  present  in  the  underlying  base  model 
described  in  the  previous  chapter.  Variable  elements 
specifying  particular  models  are  introduced  either  by 
parameter  definition  or  by  selection  of  optional  modules  in 
the  "model  pool". 

5.1  Choice  of  Programming  Environment 

Two  overriding  considerations  influenced  the  choice  of 
programming  environment  for  the  project.  Firstly,  a 
developing  system  is  a  dynamic  structure,  with  a  spectrum  of 


\ 
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activities  in  the  system  which  can  change  with  time,  and 
proliferating  cells.  Of  paramount  concern,  therefore,  is 
the  need  for  convenient  (and,  hopefully  efficient)  handling 
of  dynamically  changing  data  structures.  Secondly  the 
parallel  nature  of  the  model,  with  different  processes  going 
on  simultaneously  in  different  cells,  influences  a  decision 
concerning  possible  use  of  available  simulation  packages. 

It  has  been  said  that  the  appropriate  time  base  for  the 
base  model  is  continuous,  because  although  discrete 
developmental  events  occur,  they  arise  through  cumulative 
conditions  in  individual  cells.  Such  events  cannot 
therefore  be  scheduled  in  the  continuous  time  base  but  must 
be  recognized  as  they  occur.  The  model  is,  however,  of  the 
cell  space  type,  and  efficient  handling  of  the 
representat i on  of  parallel  processes  by  sequential  means  has 
been  accomplished  in  some  such  models  by  use  of  discrete 
event  strategies  (Zeigler,  1976;  Hogeweq ,  1978).  The 

technique  involves  processing  only  those  cells  which  have  a 
possibility  of  undergoing  change  at  the  next  time  step. 
Thus  a  "next  event"  list  is  maintained  and  languages  such  as 
SIMULA  become  useful  in  minimising  processing.  SIMULA 
moreover  simulates  some  aspects  of  parallel  processing  for 
the  user . 

However  the  benefits  of  a  discrete  event  approach  to 
cell  space  modelling  depend  on  the  level  of  activity  in  the 
cells.  It  becomes  efficient  only  if  a  significant  number  of 
cells  does  not  change  at  each  time  step.  The  developmental 
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model  has  a  high  level  of  activity;  it  is  probable  that  all 
cells  are  engaged  in  some  relevant  activity  all  the  time. 

While  packages  for  continuous  time  simulations  are 
available,  none  of  these  "differential  equation"  packages  is 
orientated  towards  cell  space  models,  and  the  type  of  data- 
handling  required  militates  against  their  use.  Typically, 
activity  in  a  continuous  system  is  modelled  by  one 
differential  equation,  and  structures  in  which  activities 
are  locally  mediated  cannot  be  treated  by  the  package. 

The  language  chosen  for  implementation  of  the  lumped 
model  is  ALGOLW.  It  has  the  right  kind  of  dynamic  data 
acquisition  handling.  In  ALGOLW  there  are  "pointer"  data 
types,  in  which  data  is  stored  and  referenced  automatically 
by  means  of  pointers,  thus  allowing  linked  allocation  of 
storage.  Moreover  the  data  thus  referenced  can  be  of  mixed 
type,  thus  allowing  a  variety  of  pertinent  information  to  be 
maintained  for  the  structure.  Pointers,  and  therefore  the 
data  structures  to  which  they  refer,  can  be  allocated  (and 
destroyed)  dynamically,  and  thus  there  is  natural  expression 
for  structural  events  such  as  cell  division.  ALGOLW  has 
efficient  numerical  capabilities,  with  an  interface  to 
FORTRAN  facilities  if  desired.  A  modelling  system 
implemented  in  this  language  should  be  widely  understandable 
and  have  adequate  flexibility. 
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5.2  Treatment  of  the  Time  Base 

The  continuous  time  base  of  the  base  model  is 
approximated  in  the  lumped  models  by  a  discrete  time  base. 
The  size  of  the  time  step  is  a  parameter  determined  by  the 
particular  model  being  studied.  For  a  given  model,  the  time 
step  may  also  be  adjusted  in  response  to  the  generated  data. 
Within  one  time  step  various  activities  may  be  taking  place, 
for  example,  cell  divisions  and  changes  in  several  state 
variable  values.  In  general  it  will  be  most  efficient  to 
complete  all  activities  for  each  cell  before  proceeding  to 
the  next  cell,  as  a  minimum  number  of  storage  operations  are 
thus  performed,  but  the  option  exists  for  processing  all 
cells  for  each  activity  in  turn  if  such  should  appear 
necessary.  Further  discussion  of  issues  concerning  the 
treatment  of  time  will  be  discussed  in  the  context  of 
particular  models  in  Chapter  Six. 

5.3  Data  Structures 

The  structural  elements  introduced  in  the  previous 
chapter  (the  individual  cell,  the  code  group,  and  the 
organism)  are  described  by  a  set  of  data  structures  in 
ALGOLW. 

5.3.1  Cells 

A  particular  cell's  internal  character  is  defined  by 
the  values  of  its  state  variables  (parameters  such  as  clonal 
markers,  local  concentration  values,  and  the  genetic  code). 
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In  the  ALGOLW  implementation  of  a  particular  model,  the 
appropriate  binary  code  string  and  parameter  fields  are  held 
for  each  cell  in  a  "state"  record  of  a  declared  structure. 
Parameters  may  include  local  concentration  of  relevant 
substances.  A  state  record  for  a  cell  is  pictured  in  Figure 
5.1.  As  cells  are  created  initially  or  through  cell 
division  such  records  can  be  allocated  through  invocation  of 
the  structure  declaration. 

When  considered  as  a  member  of  a  population,  a  cell  has 
position  and  neighbourhood  relationships.  A  two-dimensional 
population  of  cells  is  approximated  by  a  hexagonal  grid; 
that  is,  cells  are  pictured  as  being  centred  on  points 
arranged  hexagonal ly  on  a  plane.  This  choice  was  made  for 
the  following  reasons.  In  cells  under  no  distorting  forces 
on  a  plane,  spherical  closest  packing  most  nearly 
approximates  nature,  in  that  cells  appear  in  reality  to  be 
in  physical  contact  with  between  5  and  6  neighbours 
(Lawrence,  1975);  there  is  no  ambiguity  regarding 
neighbourhood  relationships  due  to  local  symmetry;  and 
arbitrary  curves  joining  cells,  as  for  example  along 
population  boundaries,  are  more  easily  approximated  than  on 
a  conventional  rectangular  array. 

The  hexagonal  grid  of  cells  is  defined  by  a  set  of 
"position"  records,  each  containing  the  hexagonal  address, 
six  pointers  to  similar  records  for  the  six  adjacent  cells, 
and  certain  other  fields.  A  position  record  for  a  cell  is 
pictured  in  Figure  5.2.  Neighbourhood  relationships  between 
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FIGURE  5.1  State  record  for  a  cell 
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cells  are  pictured  in  Figure  5.3.  One  of  these  fields  is  a 
pointer  to  a  "state"  record,  symbolizing  the  location  of  a 
particular  cell  of  individual  character  or  "name"  at  that 
hexagonal  address.  It  may  be  noted  that  if  a  cell  moves, 
for  example  in  response  to  cell  division  in  its  vicinity, 
this  movement  can  be  accomplished  by  letting  the  "position" 
record  for  its  new  position  point  to  its  "state"  record  or 
name;  thus  neighbour  links  do  not  have  to  be  changed  as  it 
assumes  a  new  position,  as  the  neighbour  links  are 
associated  with  positions  and  not  cell  names.  Other  fields 
in  the  position  record  will  be  pointers  associating  the 
position  (and,  via  the  pointer  to  the  "state"  record,  a 
particular  cell)  with  a  processing  order  and  with  a  code 
group  of  genetically  identical  cells. 

It  may  be  noted  that  in  the  formal  structure  for  a 
cell,  <T,X,W,Q,Y,f ,g>, 

the  elements  W,Q,  and  f  have  been  translated  into  computer 
terms.  W,  the  set  of  input  segments  accessing  each  cell,  is 
specified  by  the  neighbour  links  to  the  cell  and  by  a  given 
starting  condition  in  an  array.  0  is  the  set  of  state 
variables  held  in  the  cell's  "state"  record;  and  f  is 
defined  by  the  genetic  code's  switching  pattern  on  the  state 
variables.  Furthermore  structures  exist  for  selecting  Y  and 
g  in  a  given  experiment;  for  example,  the  output  function  g 
may  select  of  cells  of  identical  genetic  nature,  in  which 
case  the  list  of  cells  linked  together  by  the  group  link  can 
be  output  in  linear  or  graphic  form. 
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FIGURE  5.3  Neighbourhood  relationships  between  cells 


Positions  in  the  hexagonal  array  are  organised  into 
processing  lists  and  code  group  lists.  In  order  to  locate  a 
cell  by  its  position,  unless  further  information  is  given, 
no  alternative  to  searching  along  one  of  the  lists  exists. 
Such  searches  are  expensive  because  of  the  number  of  storage 
references  involved  in  following  the  pointers.  Whenever 
this  Kind  of  problem  arises,  local  "hash  tables"  based  on 
position  coordinates  can  be  set  up;  that  is,  cells  in  a  list 
can  be  located  in  a  linear  array  by  a  function  operating  on 
their  addresses.  The  entry  in  each  position  of  the  hash 
table  is  a  pointer  to  the  actual  cell  located  there  (the 
machine  address  of  the  record  associated  with  the  cell). 
When  no  longer  needed  these  tables  are  destroyed. 

5.3.2  Groups  of  Cells  of  Identical  Genetic  Nature 

The  formal  specification  of  a  code  group  has  been  given 
as  <T , Q , Y , f ,g> . 

The  state  variable  set  consists  of  the  genetic  code 
common  to  cells  in  the  group,  a  list  of  member  cells,  and  a 
shape  description.  The  list  of  member  cells  is  given  by 
following  the  "group  links"  in  the  position  records .  This 
list  may  therefore  be  specified  by  a  pointer  to  the  first 
cell  in  the  linked  list.  The  genetic  code  of  the  group  is 
carried  by  each  member  cell.  A  shape  description  such  as  a 
perimeter  list  can  be  constructed  for  a  given  modelling 
experiment.  For  example,  a  "shape"  record  class  may  take 
the  form  of  a  pointer  to  the  "state"  record  of  the  perimeter 
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cell  and  a  pointer  to  the  next  "shape"  record.  Then  the 
shape  list  is  given  by  a  pointer  to  its  first  record.  A 
"shape"  record  list  is  pictured  in  Figure  5.4. 

The  transition  function  f  operating  on  the  state  of  a 
code  group  must  be  specified  for  individual  developmental 
models  in  the  form  of  a  routine  which  implements  a  growth 
model  or  recognizes  critical  configurations  in  the  member 
cells.  The  output  set  Y  and  function  g  must  also  be 
specified  for  particular  requirements,  but  advantage  can  be 
taken  of  the  existing  structures  (for  example  if  Y  is  the 
perimeter  list,  g  merely  lists  it  by  following  links). 

5.3.3  The  Developmental  System 

The  developmental  system  has  been  specified  formally  by 
<T,Q,f>,  Q  being  a  list  of  code  groups  with  associated 
location  parameters.  The  list  is  defined  by  a  record  class 
known  as  "groups"  which  gives  each  group  an  identifying 
number  and  associates  with  the  number  pointers  to  the  member 
list  and  shape  description  and  appropriate  location 
parameters.  The  transition  function  must  again  be  specified 
for  individual  developmental  models.  The  "groups"  record 
list  for  a  developmental  system  is  pictured  in  Figure  5.5. 

5.4  Sequential  and  Parallel  Processing 

The  base  model  is  parallel  in  two  senses:  different 
cells  are  undergoing  changes  simultaneously,  and  different 
processes  are  occurring  at  the  same  time.  Both  kinds  of 
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parallelism  must  be  approximated  by  sequential  processing  in 
the  digital  computer. 

Parallel  processes  may  easily  be  approximated  in  a 
sequential  manner  (that  is,  one  process  completed  for  a  time 
step,  then  another)  as  long  as  they  are  independent. 
Consideration  must  be  given,  however,  to  situations  in  which 
processes  are  not  independent.  For  example  two  substances  X 
and  Y  may  be  produced  according  to  a  coupled  pair  of 
production  functions  as 'they  catalyze  or  inhibit  each  other. 
Then  the  action  of  the  transition  function  on  the  state 
variables  symbolizing  X  and  Y  must  be  implemented  so  as  to 
reflect  the  interdependence  of  their  action. 

Simultaneous  activity  in  the  array  of  cells  must  be 
approximated  sequentially  with  caution.  For  each  type  of 
activity,  care  must  be  taken  that  systematic  distortions  are 
not  introduced,  for  example  by  the  order  in  which  cells  are 
processed.  Models  of  a  particular  activity  must  be  assessed 
with  respect  to  their  success  in  preventing  such  distortion. 
In  addition  extra  storage  will  be  needed  in  order  that  the 
"simultaneity"  be  preserved  between  communicating  cells.  If 
the  next  state  of  cell  B  depends  on  the  present  state  of 
cell  A,  cell  A' s  computed  new  state  must  be  stored  until 
cell  B;  s  new  state  is  computed,  and  then  cell  A  can  be 
updated . 

Problems  of  parallel  to  sequential  conversion  will  be 
considered  in  more  detail  in  the  discussion  of  individual 
models  in  Chapter  Six. 


' 

f  .  I  I 


115 


5.5  Representation  of  Growth 

As  has  been  noted,  growth  is  an  activity  common  to  all 
developmental  models;  in  addition,  cons iderat ion  of  the 
implementation  of  growth  models  provides  an  example  of 
processing  using  the  described  data  structures.  The 
mechanics  of  cell  division  in  these  data  structures  are  is 
pictured  in  Figure  5.6. 

Cell  division  and  the  control  of  shape  in  populations 
of  cells  has  been  modelled  in  terms  of  two  decisions,  at 
each  of  which  various  model  options  can  be  brought  into 
play.  Firstly  a  cell  decides  whether  it  is  going  to  divide 
or  not  (decision  point  A).  It  does  so  according  to  some 
growth  model  (whether  an  internal  clock,  internal  threshold 
value,  external  signal,  or  probabilistic  function  value). 
If  it  does  not  divide,  processing  passes  to  another  cell. 

The  shape  of  the  population  may  be  governed  entirely  by 
this  first  decision,  for  example  if  only  those  cells  in  a 
certain  "growth  zone"  are  permitted  to  divide.  It  is 
possible,  however,  that  the  cell  may  decide  in  which 
direction  to  divide  (decision  point  B),  according  to  some 
other  model  (in  the  direction  of  least  pressure,  or  down  a 
gradient,  for  example). 

Growth  Algorithm 

(1)  Set  pointer  P  to  head  of  processing  list 

(2)  Is  P  null?  If  yes,  stop. 

(3)  Decision  point  A:  does  cell  (B)  divide?  If  yes, 

division  model.  If  no,  go  to  9. 


i nvoke 


. 

' 
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(4)  Decision  point  B:  In  what  direction  does  it  divide? 
Chose  random  direction,  or  invoke  direction  model, 
setting  direction  i,  Ki<6. 

(5)  Proceed  in  "position"  data  structure  along  neighbour 
link  i  until  first  empty  space  is  reached  (link  is 
null). 

(6)  Create  new  "position"  record  for  empty  space  and  link  to 

appropriate  neighbours. 

(7)  Move  "state"  record  pointer,  processing  links  and  group 
links  outwards  along  direction  i  to  symbolize  expansion 
of  population  making  room  for  new  cell. 

(8)  Create  new  "state"  record  for  new  cell,  and  attach  to 

adjacent  position  i.  Put  it  before  cell  (P)  on 
processing  list  (it  should  not  divide  in  this  time 
step)  and  link  it  to  the  same  code  group  as  cell  (P). 

(9)  Set  P  to  the  link  of  P  (go  to  next  cell  on  processing 
list).  Go  to  2 . 

5.6  Problems  of  Economization 

The  data  structure  necessary  to  implement  a  model  in 
which  cells  carry  internal  information,  and  in  which 
processing  is  efficient,  is  unavoidably  extensive.  It  is 
likely  that  in  an  experiment  of  some  size  the  space 
requirements  will  become  too  large.  There  are  three  ways  in 
which  economization  of  space  could  be  realized  when 
necessary . 

The  first  technique  involves  adjusting  the  focus  of 
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attention.  For  example  during  simulation  of  a  disc 
undergoing  compar tmenta 1 i zat i on ,  different  code  groups  will 
be  formed  cor respondi ng  to  the  different  compartments.  If 
the  disc  were  to  become  too  large  to  handle  in  entirety,  the 
simulation  could  be  stopped  and  restarted  with  focus  on  one 
compartment.  In  this  case  the  entire  processing  list  would 
be  replaced  by  a  code  group  list. 

The  second  technique  involves  representing  a  number  of 
cells  by  one  cell  with  representat i ve  characteristics.  A 
mapping  of  the  neighbourhood  geometry  from  the  simple  case 
to  the  averaging  case  would  have  to  be  preserved,  and  the 
method  would  not  be  applicable  in  cases  where  cells  in 
different  code  groups  were  i nter sper sed . 

The  third  technique  involves  a  redefinition  of  the  data 
structures  used  in  implementing  the  lumped  model.  As  they 
have  been  defined,  the  data  structures  mirror  the  structures 
defined  for  the  model  specification  in  Chapter  4.  However  a 
good  deal  of  space,  and  some  time  spent  in  storage 
references,  could  be  saved  by  making  greater  use  of  hash 
tables . 

Before  a  processing  step,  comprising  cell  division, 
communication  and  possible  genetic  events,  the  size  of  the 
organism  is  known  by  parameters  such  as  the  extent 
parameters  given  for  each  code  group.  An  estimate  can  be 
made  of  the  probable  size  the  organism  will  attain  during 
this  step  and  a  hash  table  set  up  to  cover  it.  This  hash 
table  would  contain,  in  each  row  cor  respond i ng  to  each 


' 


. 


. 

' 


. 


119 


position,  all  of  the  structures  as  previously  defined  except 
the  "state"  record.  It  would  contain  fields  for  the 
processing  link  (row  number  of  the  next  cell  to  be 
processed),  group  link  (row  number  of  another  cell  in  the 
same  group),  "state"  record  link  (pointer  to  the  record)  and 
perimeter  link  for  shape  description  (null  for  internal 
cell,  row  number  of  next  perimeter  cell  for  external  cell). 
Neighbourhood  links  would  be  obtainable  by  arithmetic  on  the 
row  numbers  of  the  hash  table. 

The  hash  table  representat i on  is  more  economical  but 
does  not  mirror  model  structures  as  closely  as  the  record 
classes  defined  previously.  Processing  could  be  conceivably 
less  convenient.  Until  the  modelling  framework  has 
undergone  the  first  developmental  phase,  the  previously 

defined  data  structures  will  be  maintained. 

5.7  Validity  of  the  Implementation 

Specification  morphisms  between  the  computer 
implementation  and  the  lumped  models  are  determined  by 
inspection  of  the  state  transition  function.  They  cannot  be 
discussed,  therefore,  until  specific  developmental  models 
with  specific  transition  functions  are  defined.  Care  has 
been  taken,  however,  to  ensure  a  degree  of  structural 

faithfulness  in  defining  data  structures  to  mirror  the 

structural  elements  of  the  underlying  lumped  model,  namely 

the  cell,  the  code  group,  and  the  developmental  system. 
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CHAPTER  SIX 

Experiments  and  Results 


Thi  s 

chapter  describes  the 

construct 

ion  and 

use 

of  a 

mode 

1  1  ing 

structure  containing 

many  opt 

ions  and 

strategies 

for 

simulating  development. 

In  the 

context 

of 

thi  s 

structure,  an  "experiment"  is  the  simulation  of  a  particular 
model  in  a  suitable  experimental  frame.  The  establishment 
of  the  entire  structure,  experiment  by  experiment,  can  in  no 
sense  be  an  exhaustive  process;  that  is,  there  are  always 
further  options  to  be  explored  in  modelling  a  given 
phenomenon.  Attention  is  given  to  models  of  current  genetic 
interest,  for  example,  reaction-diffusion  systems; 
strategies  for  implementing  other  models  are  mentioned  where 
re  1 evant . 


6.1  Experimental  Sequence 

As  previously  discussed,  development  involves  the 
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interaction  of  growth  (cell  divisions  in  a  population  of 
cells)  and  pattern  formation.  In  the  initial  stages  of  the 
project,  attention  was  focussed  on  experiments  with  models 
of  growth.  Subsequently  attention  shifted  to  investigation 
of  mechanisms  of  communication  between  cells  and  the 
formation  of  patterns  of  information  in  groups  of  cells. 
Experiments  were  then  aimed  at  modelling  growth  and 
intercellular  communication  simultaneously.  Finally, 
various  models  for  actual  developmental  systems  were 
studied . 

6.2  Modelling  Growth 

Growth  may  be  defined  informally  as  the  process  by 
which  a  developmental  system  is  enlarged  by  division  of 
individual  cells.  The  process  is  controlled  at  both 
cellular  and  system-wide  levels.  Each  cell  appears  to 
divide  in  response  to  attainment  of  a  critical  internal 
level  of  protein  content  (Alberghina,  1978).  Cell  division 
rate  in  a  system  appears  to  be  constant  in  a  given 
environment  and  to  be  responsive  to  changes  in  environment 
(Kiefer,  1968;  Wheldon,  1973).  Furthermore,  the  Kinetics  of 
cell  division  may  alter,  both  in  response  to  cell 
differentiation  and  in  response  to  attainment  of  maximum 
allowable  size  of  the  system  (Alberghina,  1975,  1977;  Kim  et 
al.,  1973).  It  may  be  remarked  therefore  that  models  of 
growth  should  embody  various  levels  of  control  and  feedback. 
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6.2.1  Formal  Specification  of  Growth  Models 

A  growth  model  is  a  structure  G  given  by 
<Q,T,y>,  where 

Q  is  the  single  state  variable,  consisting  of  a  linked  list 
of  cells  defined  by  their  addresses,  links  to  neighbours, 
and  pointers  to  individual  parameters; 

T  is  a  continuous  time  base; 
y  is  the  state  transition  function 
y : Q  -->  Q. 

The  state  transition  function  y,  which  maps  a  list  of 
cells  into  another  list  containing  new  daughter  cells, 
operates  in  the  continuous  time  base  T.  Its  action  is  to  be 
represented  in  discrete  time,  and  in  iterative  form,  in  a 
lumped  model  given  by  systems  SI  and  S2  at  the  cellular 
level,  and  G1  at  the  group  level. 

SI  is  a  system  <Q 1 , I , T 1 , f 1 , g 1 , Y 1 > ,  where 
Q1  is  the  state  variable  set,  each  element  of  which  consists 
of  the  address  of  a  cell; 

I  is  the  set  of  all  inputs  to  SI  (factors  in  the  cell's 
environment  affecting  division); 

T1  is  a  discrete  time  base,  with  a  constant  step  size 
measuring  the  time  to  consider  one  cell  (cell  processing 
time) ; 

fl  is  the  state  transition  function,  which  replaces  one  cell 
by  the  next  on  the  processing  list 
f 1 : Q 1  -->  Q 1  is  given  by 

f 1 ( c ) = 1 i nk ( c )  ,  c  6  Q1.  If  c=null,  system  terminates. 
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gl  is  the  output  function,  discussed  below; 

Y1  is  the  output  set,  consisting  of  all  couplets  (a,b), 
a  6  {0,1},  be  11,2,3,4,5,6}. 

The  output  function  g 1 : Q 1  x  I  -->  Y,  calculates  a 
condition  local  to  the  cell  in  question  and  specific  to  a 
particular  model.  It  returns  a=0  if  the  cell  is  not  to 
divide,  and  a=1  if  it  is.  Furthermore  it  returns  the  digit 
b  signifying  the  direction  of  cell  division,  which  is  again 
specific  to  a  particular  model. 

The  output  couplets  of  SI  form  input  to  subsystem  S2, 
which  is  a  system 

<Q2 , 1 2 , T 1 , f 2 ,g2 , Y2> . 

The  system  S2  is  thus  a  means  of  locating  dividing 
cells  on  a  processing  list, 
where 

Q2  is  the  single  state  variable  set,  with  domain  the 
positive  integers; 

12  is  the  input  set,  generated  as  Y1  for  SI; 

T1  is  the  discrete  time  base  of  system  si; 
f2  is  the  state  transition  function 

f2:Q2  x  12  -->  Q2  where  f2(q,a)=q+1; 
g2  is  the  output  function;  it  generates  a  1  if  a  dividing 
cell  is  reached  in  the  processing  list,  0  otherwise: 
g2:Q2  x  12  -->  Y2  where 
g2 (q , a ) =0 ,  a=0 
g2 ( q , a  )  =  1  ,  a= 1 ; 

Y2  is  the  set  (0,1). 
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The  system  G1  is  given  by 
<Q , I G , XG , T 1 , f G> ,  where 

Q  is  the  single  state  variable,  consisting  of  a  linked  list 
of  cells  with  addresses  and  neighbour  links; 

IG  is  the  input  set,  a  triplet  (c,d,e),  where 
c  is  the  output  of  S2, 
d  is  the  state  variable  value  Q2  of  S2, 
e  is  b  from  system  SI  (direction  of  division); 

XG  is  the  set  of  actual  values  received  from  $2 ; 

T1  is  the  discrete  time  base; 
fG  is  the  state  transition  function 
f G : Q  x  IG  -->  Q  where 
f G (q , 0 , d , e ) =q  ,  q  6  Q 
fG(q,1,d,e)  results  in  the  action 

(a)  copy  entry  d+1  in  the  list  q  and  insert  it  in  position 
d+2  in  the  list  (set  pointers  to  and  from  it); 

(b)  change  neighbour  links  where  necessary; 

(c)  change  parameter  pointers  along  direction  e  effecting 
cell  movements. 

The  operation  of  the  growth  model  may  be  summarized 
verbally  as  follows.  A  structure  of  cells,  specified  by  a 
linked  processing  list  with  neighbour  links,  undergoes 
changes  which  are  mediated  dynamically  at  the  level  of  the 
individual  cell.  At  discrete  intervals  of  time  the 
structure  changes  by  the  addition  of  a  single  cell  and 
cor  respond i ng  changes  in  position  in  the  surrounding 
structure.  Growth  in  a  given  time  interval  is  represented 
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by  the  sum  of  such  individual  events  arising  dynamically 
during  that  time. 

Implicit  in  the  formal  presentation  of  the  growth  model 
are  concepts  developed  by  Zeig ler ( 1 976 )  in  his  Twelve 
Postulates  for  valid  modelling  and  simulation.  In 
particular,  the  "real  system"  is  the  set  of  all  cellular 
structures  in  nature  undergoing  cell  divisions.  It  is 
observable  through  data  such  as  initial  and  final  shape, 
number  of  cells,  and  clonal  distributions  (marked  cell 
lineages).  The  "base  model"  is  the  structure  G  above.  The 
set  of  "experimental  frames"  in  which  G  could  operate 
includes  all  cell  structures;  however,  G  is  studied  here  in 
the  universe  of  two-dimensional  structures  or  organisms 
observable  in  two-dimensional  cross-section.  A  particular 
experimental  frame  is  given  by  a  set  of  allowable  input 
values  (a  subset  of  the  environmental  influences  of  set  I  of 
SI)  and  a  set  of  observable  outputs  (a  subset  of  the 
observables  mentioned  above  in  the  "real  system").  The 
"base  model  in  the  experimental  frame",  G/E,  is  the 
structure  G  with  Q  restricted  to  those  states  for  which  the 
observables  fall  within  the  output  set  of  the  experimental 
frame,  an  output  set  consisting  of  the  output  set  of  the 
experimental  frame  plus  a  null  output  for  nonvalid  (in  E) 
input  segments,  the  state  transition  function  f  restricted 
to  final  states  which  have  observables  in  the  output  set  of 
E,  and  an  output  function  yielding  these  observables  for 
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The  operation  of  the  base  model  in  E  (a  particular 
experiment)  is  Known  to  the  observer  through  observations 
made  on  the  real  system.  The  real  system  itself  is 
accessible  through  the  totality  of  observations  in  all 
experimental  frames. 

In  the  process  of  modelling  and  simulating  growth,  two 
iterative  specifications  are  developed.  Firstly,  there  is 
the  lumped  model  (systems  SI,  S2  and  G1),  and,  secondly,  the 
simulating  program.  There  are  two  tests  of  validity  of  the 
process  which  must  be  made.  If  the  data  generated  by  the 
lumped  model  in  experimental  frame  E  are  the  same  as  the 
data  expected  for  the  base  model  in  E,  then  the  lumped  model 
is  said  to  be  valid  for  the  real  system  in  E  (the  action  of 
SI,  S2  and  G1  on  allowable  environmental  inputs  and  initial 
structures  produces  realistic  final  structures  for  this 
experimental  frame).  If  there  is  a  specification  morphism 
between  the  simulating  program  and  the  lumped  model,  then 
the  program  is  said  to  be  a  valid  simulator  of  the  lumped 
model  . 

The  success  of  the  iterative  lumped  model  in  portraying 
reality  can  only  be  judged,  therefore,  by  comparing  the 
iteratively  generated  results  with  real  data.  The  success 
of  the  simulating  program  in  translating  the  intentions  of 
the  lumped  model  can  be  ensured  by  requiring  that  (1)  the 
input  segments  of  the  computer  program  can  all  be 
"translated"  into  allowable  input  segments  of  the  lumped 
model  (thus  ensuring  that  extra  influences  are  not  at  work 
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in  the  simulator),  (2)  there  is  an  "onto"  map  from  a  subset 
of  the  state  set  of  the  lumped  model  to  the  state  set  of  the 
simulator  (every  state  of  the  simulator  represents  a  valid 
state  in  the  lumped  model)  and  (3)  any  output  of  the  lumped 
model  can  be  identified  with  some  output  of  the  simulator 
(ensuring  that  the  simulator  has  a  broad  enough  scope). 
Then  if  both  the  state  transitions  and  outputs  are  preserved 
under  translation,  a  specification  morphism  holds. 

In  terms  of  the  modelling  of  growth,  the  computer 
program  is  designed,  for  each  of  the  options  developed 
below,  so  as  to  accept  allowable  inputs  and  initial  states 
for  the  given  experimental  frame.  The  actual  state  values 
which  may  be  generated  by  the  program  (structures  of  cells) 
are  restricted  to  structures  realizable  in  the  lumped  model 
(for  example,  by  restricting  the  structures  to  the  plane). 
Requirement  (3),  that  any  output  of  the  lumped  model  can  be 
achieved  by  the  simulator,  is  harder  to  ensure.  It  can  be 
inferred  only  after  a  multiplicity  of  simulations  has  been 
run.  Preservation  of  the  state  transition  function  means 
that  when  the  representat ion  of  a  state  of  the  lumped  model 
in  the  simulator  undergoes  transition,  the  resultant  state 
of  the  simulator  represents  that  particular  state  of  the 
lumped  model  under  the  same  transition;  that  is,  a  structure 
generated  by  the  simulator  representat i on  of  the  transition 
function  itself  represents  that  structure  which  would  have 
been  generated  by  the  transition  function  of  the  lumped 
model.  That  this  condition  holds  must  be  established  by 
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simulation  experiments  (see  Sections  6. 2. 3-6. 2. 5  below). 
Similarly,  output  preservation  is  established  by  running 
s i mu  1  at i ons . 

The  major  concern  in  establishing  the  validity  of  the 
lumped  (iterative)  model  as  a  representat i on  of  the  real 
system  lies  in  the  replacement  of  a  continuous  time  base  by 
a  discrete  one.  Thus  Section  6.2.2  is  concerned  with 
experiments  designed  to  validate  this  representation. 

Before  particular  experiments  in  growth  are  described, 
some  general  introductory  remarks  may  be  made.  Firstly,  the 
model  of  growth  is  necessarily  complex  and  mul t i level  led 
because  it  deals  with  dynamically  changing  structure,  one  of 
the  essential  challenges  in  the  project,  and  because  it  must 
contain  opportunities  for  different  levels  of  control. 
Secondly,  the  system  G  is  represented  as  being  autonomous 
( i nput - f ree ) .  Environmental  influences  are  modelled  at  the 
level  of  the  individual  cell  in  system  SI,  reflecting  the 
fundamental  assumption  that  the  processes  of  development  are 
mediated  by  events  local  to  cells.  Thirdly,  it  may  be  noted 
that  several  different  time  bases  are  employed  in  the  model. 
The  relationship  between  these  time  bases  should  be 
explored . 

6.2.2  Sequential  and  Parallel  Processing:  Experiments  with 
T  i  me 

The  system  G  proceeds  from  one  value  of  its  time  base  T 
to  the  next,  producing  a  new  structure  of  cells  and  daughter 
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cells.  In  nature,  individual  cell  divisions,  once  initiated 
appear  to  proceed  independently  without  regard  to  other 
divisions.  That  is,  in  the  context  of  a  single  interval  of 
time  of  base  T,  all  cells  that  divide  do  so  simultaneously. 
System  G  could  therefore  be  modelled  with  structural 
validity  by  means  of  parallel  processing.  However  systems 
SI,  S2  and  G1,  which  are  iterative  specifications  suited  to 
a  digital  computer,  treat  the  cells  sequentially  (as  they 
appear  on  a  processing  list).  An  hypothesis  concerning  the 
validity  of  such  processing  is  presented  and  tested  at 
various  levels. 

Processing  Hypothesis:  The  system  G  may  be  validly  modelled 
by  sequentially  processing  SI,  followed  by  S2  and  G1. 

The  processing  hypothesis  may  be  tested  at  two  levels. 
Firstly,  the  effects  of  separating  the  decision  to  divide 
(SI)  from  division  itself  (S2  and  G1)  should  be 
investigated.  Secondly,  any  distortions  in  structure 
resulting  from  cells  dividing  in  sequence  (in  S2  and  G1) 
should  be  noted. 

Separation  of  SI  from  (S2  and  G1)  is  desirable 
computationally  because  cells  already  processed  in  S2  and  G1 
are  thus  incapacitated  from  dividing  in  this  time  step  in 
the  base  T1.  If  they  were  allowed  to  decide  to  divide  while 
the  division  process  is  in  train,  the  list  would  have  to  be 
re-scanned  to  accomplish  their  division,  and  the  process 
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would  repeat  itself  indefinitely.  The  question  is  whether 
this  limitation  of  the  model  constitutes  a  serious 
difference  from  the  natural  situation  being  modelled.  In 
nature,  during  an  interval  of  time  (which  is  continuous), 
cells  initiate  the  division  process  at  any  time. 

Thereafter,  however,  until  the  process  is  complete  division 
cannot  again  be  initiated.  Thus  if  the  discrete  time  step 
in  system  G  is  not  allowed  to  be  greater  than  the  cell  cycle 
time,  no  cells  "flagged"  to  divide  in  that  time  step  will  do 
so  again.  However  cells  not  initially  flagged  to  divide  may 
subsequently  decide  to  do  so  because  of  a  change  in  local 
conditions  (modelled  in  SI  by  the  action  of  gl  in  response 
to  the  state  of  the  cell  itself  and  inputs  from  its 
environment.)  If  any  discrete  system  modelling  such  change 
in  local  conditions,  is  not  allowed  to  have  a  time  step 
incremented  during  a  time  step  of  G  (that  is,  the  time  step 
of  that  system  is  an  integral  multiple  of  the  time  step  of 
G),  such  a  change  cannot  happen. 

Thus  there  are  conditions  on  the  time  base  of  G  which 
allow  it  to  be  modelled  by  separation  of  SI  from  S2  and  Gl. 
In  terms  of  the  relationship  of  the  lumped  model  of  G  to  a 
real  system  it  models,  the  conditions  under  which  a  discrete 
time  base  can  represent  continuous  time  are  beginning  to  be 
established:  the  time  step  must  be  smaller  than  cell  cycle 

time,  and  "small  enough"  to  be  sensitive  to  other  processes 
influencing  cell  growth. 

The  above  discussion  of  separating  sequential 
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processing  into  SI  and  (S2  and  G1 )  has  introduced  the 
question  of  changes  in  local  conditions  which  are  effected 
by  systems  other  than  the  growth  model.  The  question  cannot 
therefore  be  pursued  experimentally  without  the  development 
of  such  systems.  Restrictions  on  the  time  base  have  been 
established,  and  attention  is  now  turned  to  the  second  level 
at  which  the  processing  hypothesis  can  be  tested.  An 
experiment  designed  to  test  any  distortion  of  structure  by 
sequential  division  of  cells  is  presented. 

Processing  Experiment:  To  examine  distortions  of  structure 
arising  from  sequential  processing. 


Experimental  frame:  Planar  structures  are  constructed  using 
the  growth  model  in  various  forms.  The  initial  structures 
are  identical  but  processing  is  undertaken  in  varied 
ordering.  Autonomous  models  are  used  (growth  unaffected  by 
environment,  no  input),  and  output  is  to  consist  of  a 
mapping  of  the  structures  at  each  growth  time  step,  where 
cells  bear  "clonal"  markers  -  that  is,  the  lineage  of  each 
cell  is  given . 


Model  features:  The  output  function  gl  of  SI 


are  flagged  to  divide) 
synchronous  division), 
direction,  is  "uniform 
number  generator)  or 


g i ves  1  for  all  cells 
Its  second  component, 
random"  (generated  by 
in  the  direction  of  " 


(by  which  cells 
( uncond i t i ona 1 , 
the  division 
a  pseudo- random 
least  pressure" 


* 

b  *«  :rr  i  ■ 


132 


(generated  by  looking  along  the  six  neighbouring  directions 
and  finding  the  nearest  empty  space) . 

The  transition  function  f2  of  S2 ,  which  follows  the 
processing  list  and  replaces  one  cell  by  the  next,  is  varied 
by  means  of  manipulating  the  processing  list.  This  list  can 
be,  in  turn,  left  in  its  natural  order  as  daughter  cells  are 
added,  or  "randomized"  at  each  growth  time  step,  or  sorted 
consistently  into  the  same  geometr i ca 1 ly-or i ented  order. 

Test  Cases:  Initial  configurations  of  various  types  are 
tested.  In  particular,  a  structure  which,  under  the  "least 
pressure"  option,  may  be  expected  to  develop  into  a 

structure  with  coherent  clones,  is  chosen;  it  consists  of 
six  cells  surrounding  a  seventh  in  hexagonal  symmetry. 
Various  irregular  structures  are  studied;  under  the  random 
direction  option,  distortion  might  be  suspected  if  cells  of 
the  same  lineage  form  unnaturally  large  coherent  clumps.  In 
all,  about  20  structures  thought  to  be  of  interest  are  taken 
through  about  8  division  cycles. 

Results:  The  output  from  the  hexagonal ly  symmetric  initial 
structure  is  presented  in  Figures  6. 1-6.3,  as  it  is 
representat i ve  of  all  results  obtained.  In  no  case  are 
completely  coherent  clones  obtained,  but  there  is  no 
striking,  systematic  difference  between  results  obtained 
with  unsorted,  randomized,  or  geometrically  ordered 
processing.  (It  should  be  noted  that  complete  coherence 
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FIGURE  6.1  Growth  with  randomized  processing  list 
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FIGURE  6.2  Growth  with  unsorted  processing  list 
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FIGURE  6.3  Growth  with  geometrically  ordered  processing  list 
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would  not  be  expected  even  with  parallel  processing  in  a 
hexagonally  symmetric  situation  as  there  is  symmetry  in  the 
choice  of  division  direction  along  the  axes  of  symmetry). 
Most  strikingly,  persistent  geometric  ordering  does  not 
appear  to  result  in  a  distortion  in  this  case.  It  might  be 
expected  that  in  Figure  6.3,  where  processing  always 
proceeds  from  the  bottom  left,  the  cell  lineages  5  and  6 
would  be  more  fragmented  due  to  cell  movements  caused  by 
later  cell  divisions;  such  is  not  the  case. 

Whatever  the  processing  order,  cells  dividing  towards 
the  direction  of  least  pressure  tend  to  shift  irregular 
initial  configurations  towards  similar  convex  ones  (Figure. 
6.4).  In  order  to  emphasize  this  result,  the  perimeter  of 
the  structure  is  displayed. 

Conclusions:  The  processing  order  in  the  sequential 

modelling  of  growth  does  not  appear  to  distort  resulting 
structures.  (The  option  of  randomizing  the  processing  list 
is  retained  in  case  systematic  distortion  is  ever 
suspected).  Based  on  the  results  of  the  processing 
experiment,  and  under  the  above  restrictions  on  the  time 
base,  the  processing  hypothesis  is  considered  established. 

Two  conjectures  may  be  advanced  on  the  strengths  of  the 
results  of  the  processing  experiment.  Firstly,  it  may  be 
that  distortions  do  not  appear  because  they  are  "damped"  by 
other  influences  which  model  "dissipative"  forces  (for 
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FIGURE  6.4  Growth  from  irregularity  to  convexity 
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example,  division  in  the  direction  of  least  pressure,  which 
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structure  implies  expending  energy.  (Care  must  be  taken  if 
such  disruptive  influences  are  modelled,  that  distortions 
from  sequential  processing  do  not  then  appear).  Secondly, 
it  may  be  suggested  that  in  modelling  systems  in  which 
coherent  clones  appear,  processes  preserving  coherence 
beyond  the  state  of  initial  contiguity  must  be  included. 
These  might  include  imposing  a  "bond"  between  cells  in  the 
same  clone  lessening  the  likelihood  of  their  being  separated 
by  other  cell  divisions.  Such  bonds  could  be  modelled  by 
recognizing  neighbouring  cells  of  the  same  clone  and 
preferring  cell  movements  which  do  not  separate  them. 

6.2.3  Growth  of  Competing  Structures 

When  the  scope  of  a  modelling  experiment  is  such  that 
different  groups  of  cells  are  processed  separately, 
conflicts  may  arise  between  groups  of  cells  growing  in 
physical  proximity.  An  example  is  compartments  growing  in  a 
Drosophila  disc.  Various  options  have  been  implemented  for 
dealing  with  structures  competing  for  space  to  expand. 

The  conflict  is  recognized  when  a  cell  selects  a 
division  direction  and  attempts  to  move  cells  outwards  along 
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that  direction.  In  searching  along  the  "ray"  for  an  empty 
space,  it  finds  instead  a  cell  belonging  to  another  group. 
Subsequently  one  of  the  following  options  may  be  selected: 
Option  1:  The  second  group  is  regarded  as  penetrable  and 
cell  movement  outward  along  the  ray  is  continued  to  the 
opposite  border  of  the  second  group. 

Option  2:  The  second  group  is  regarded  as  impenetrable  and 
as  exerting  a  force  opposed  to  growth  in  that  direction.  In 
this  case,  the  direction  of  cell  division  Is  reversed. 

Option  3:  The  second  group  is  regarded  as  impenetrable,  but 
the  cell  divides  towards  it  causing  a  movement  of  cells 
sideways  along  the  boundary  to  accomodate  the  new  cell. 

These  options  have  been  implemented  but  not  extensively 
tested,  as  they  do  not  come  within  the  mainstream  of 
development  of  the  project.  Figures  6.5  and  6.6  illustrate 
the  effect  of  a  barrier,  either  of  other  cells  or  of  a 
physical  nature,  on  a  dividing  population  under  options  2 
and  3  respectively. 

6.2.4  Modelling  the  Division  Decision 

The  output  function  gl  of  SI  generates  a  1  for  each 
cell  that  is  to  divide.  It  does  so  by  some  calculation  on 
local  conditions  (the  values  of  state  parameters  and 
inputs).  The  nature  of  this  calculation  is  dependent  on  the 
particular  model  of  growth  selected.  Autonomous, 

unconditional  division  has  been  presented  in  Experiment  A. 
Other  options  are  presented  here.  In  what  follows  "random" 
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FIGURE  6.5  Growth  near  an  impenetrable  barrier 
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means  "uniform  random"  unless  otherwise  specified. 

Option  1:  Cell  division  is  a  determinate  event  occurring  at 
regular  intervals  for  each  cell.  In  the  model,  gl  outputs  a 
1  whenever  an  internal  clock  of  the  cell  (one  of  its  state 
variables)  reaches  a  maximum  value  (measuring  the  end  of  the 
cel  1  cycle) . 

Output  from  a  particular  experiment  using  this  option 
is  presented  in  Figure  6.7.  An  initially  "random"  setting 
of  clock  values  (positions  in  the  cell  cycle)  was  imposed  on 
the  cells.  Cells  attaining  a  clock  value  of  7,  in  this 
example,  will  divide  in  the  next  growth  cycle,  and  their 
clocks  will  be  reset  to  1. 

Option  2:  Cells  divide  in  response  to  the  attainment  of  a 
threshold  value  of  a  local  function  (one  of  the  state 
variables).  An  experiment  using  this  option,  in  which  the 
part  of  the  function  is  played  by  f=x,  where  x  is  one  of  the 
position  coordinates,  with  threshold  value  of  2,  is 
presented  in  Figure  6.8  (x  is  the  horizontal  coordinate). 
It  should  be  noted  that  this  "function",  which  implies  a 
ce 1 T  s  knowledge  of  its  coordinates  (non-local  knowledge), 
is  not  intended  to  be  realistic,  but  is  merely  convenient 
for  purposes  of  demonstration. 

Option  3:  The  control  of  cell  division  has  deterministic  and 
probabilistic  elements.  Once  the  events  leading  up  to  cell 
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FIGURE  6.7  Successive  time  steps  in  clocked  growth 
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division  (protein  synthesis  and  DNA  replication)  are 
initiated,  division  will  occur  after  a  determined  interval. 


The  initiation  of 

these 

events  is 

apparently  random,  in 

propor  t i on 

to  a 

rate 

i nf 1 uenced 

by  the  environment,  and 

aborted  in 

cells 

which 

are  i  n 

the  "refractory  period" 

fol lowing 

recent 

ce  1  1 

division. 

This  model  of  growth  in 

cell  populations  is  discussed  in  Alberghina  (1978)  and  is 
similar  to  many  others  (Kim,  1975;  Kiefer,  1968;  Wheldon, 
1973) . 

Those  intracellular  controls  on  cell  division  rate 
which  are  affected  by  the  environment  are  modelled  here,  in 
an  i nterce 1 1 u 1 ar  context,  by  a  single  parameter,  the 
doubling  time  (the  time  in  which  the  population  doubles  in 
number).  The  doubling  time  is  given  as  an  integral  number 
of  time  steps  in  this  context.  In  ideal  conditions  the 
doubling  time  will  be  low,  whereas  in  certain  conditions 
(for  example,  lack  of  nutrient  species)  it  will  approach 
infinity  (a  "resting  state"). 

The  output  function  gl  of  SI  operates  here  as  follows: 

1.  Determine  what  proportion  of  cells  divide  in  this  time 
step  (t):  The  number  of  dividing  cells  ( nd )  is  (n*t/dt), 
where  n  is  the  number  of  cells  and  dt  the  doubling  time. 

2.  Choose  nd  cells  at  "random"  from  the  processing  list. 
Compare  their  internal  clocks  to  see  if  they  are  still  in  a 
refractory  state  following  prior  division;  if  so,  replace 
them  ( at  random) . 

3.  Let  gl  generate  a  1  for  these  cells. 
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Figure  6.9  illustrates  growth  using  the  random 
selection,  doubling  time,  refractory  period  model.  It 
should  be  noted  that,  for  a  doubling  time  of  2  time  steps, 
the  number  of  cells  generated  by  the  algorithm  is  within  5% 
of  the  expected  number  (n*2**k,  where  k  is  the  number  of 
times  the  doubling  period  has  been  repeated).  The 
refractory  period  is  two  time  steps  in  this  example.  In 
terms  of  an  input-output  observation  morphism,  the  model  is 
a  valid  lumped  model  of  the  Alberghina  intracellular  model. 

The  three  options  of  clocked  growth,  division  in 
response  to  a  local  function  and  the  Alberghina  model  are 
available  and  illustrate  the  division  of  individual  cells  in 
response  to  local  conditions.  Control  mechanisms  which  have 
been  modelled  include  local  ones  (internal  cell  cycle  clock) 
and  environmental  input  (doubling  time).  However,  some  of 
the  central  issues  in  controlling  growth  are  beyond  the 
scope  of  such  controls.  The  control  of  shape  (in  formation 
of  a  limb),  for  example,  and  the  limitation  of  cell  division 
once  a  critical  size  has  been  reached,  are  phenomena  which 
appear  to  have  their  domain  in  groups  of  cells  and  are 
therefore  not  simply  explicable  in  local  terms.  Under  the 
assumption  that  such  phenomena  appear  to  the  individual  cell 
in  terms  of  local  controls,  the  "global"  control  of  shape 
and  size  is  most  probably  effected  by  configurations  arising 
in  groups  of  cells  due  to  i nterce 1 1 u 1 ar  processes  or 
" commun i cat i on " . 

The  apparently  random  initiation  of  cell  division  by 
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FIGURE  6.9b  10  growth  cycles  from  7  initial  cells 
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cells  not  in  the  refractory  period  may  hide  a  varying  signal 
from  an  initiator  substance,  whose  recognition  is 
complicated  by  the  internal  clock  states  of  the  target 
cells.  Such  oscillation  phenomena  have  been  noted  in 
organisms  (Goodwin,  1973,  Reiner,  1968).  At  the  critical 
size  such  a  signal  might  no  longer  exceed  a  threshold 
initiating  cell  division.  Reaction-diffusion  systems  exist 
which  generate  travelling  waves  in  response  to  an 
oscillating  source;  for  particular  values  of  the  system 
parameters,  such  waves  are  transmitted  unattenuated  (Gmitro 
and  Scriven,  1965).  The  "critical  size"  might  be  one  in 
which  unattenuated  transmission  is  no  longer  supported. 

In  some  instances  capacity  for  cell  division  appears  to 
vanish  as  cells  become  differentiated;  however,  the  healing 
of  wounds  occurring  in  structures  of  such  cells  (Bryant, 
1971)  might  indicate  an  ability  to  respond  to  disruption  of 
a  stable  configuration  inhibiting  cell  division.  The 
differentiated  cells  either  sense  such  disruption  (perhaps 
the  critical  size  is  no  longer  operative)  with  renewed 
activity  or  are  somehow  re-enabled  to  divide.  The  latter 
alternative,  implying  an  undoing  of  the  differentiation 
process,  appears  unlikely. 

Such  global  mechanisms  of  controlling  growth  should  be 
susceptible  to  modelling  in  terms  of  division  in  response  to 
a  local  threshold  (Option  2  above),  where  the  concentration 
of  initiator  substance  is  modelled  in  terms  of  intercellular 
communication  of  the  reaction-diffusion  type. 


' 


' 


150 


6.2.5  Modelling  the  Division  Direction 

The  shape  of  a  growing  group  of  cells  clearly  can  be 
affected  by  the  division  direction  of  cells.  The  options 
mentioned  in  the  processing  experiment  were  a  random  choice 
of  direction,  division  towards  the  direction  of  least 
pressure,  and  division  down  a  concentration  gradient.  The 
problem  is  that  of  the  orientation  of  a  linearly  distributed 
mass  in  the  presence  or  absence  of  body  forces.  A  lone  cell 
might  be  expected  to  divide  in  a  random  direction  in  the 
absence  of  other  influences;  a  cell  in  a  mass  of  cells  might 
be  expected  to  experience  pressure  and  divide  in  the 
direction  of  least  pressure,  unless  more  dramatic  influences 
were  present.  However,  if  there  is  pronounced  flow  of  a 
chemical  species  in  a  particular  direction,  the  division 
axis  might  be  expected  to  lie  along  this  direction  in  the 
manner  of  a  log  floating  in  a  current.  Such  an  event  may 
occur  when  a  body  of  cells  is  entering  or  leaving  a  critical 
size  for  some  reaction-diffusion  system;  in  the  transition 
from  one  stable  state  to  another,  flow  of  the  relevant 
species  will  occur,  and  it  is  unlikely  that  any  other 
systems  have  the  same  critical  sizes  and  are  in  transition 
simul taneous ly . 

This  mechanism  suggests  a  purely  mechanical  coupling 
between  i nterce 1 1 u 1 ar  communication  and  the  control  of  shape 
in  growth.  In  a  structure  such  as  a  limb,  at  non-critical 
phases  general  enlargement  might  be  expected,  while  at 
critical  values  of  limb  length,  when  flow  up  or  down  the 
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limb  occurs,  an  elongation  phase  might  occur.  Such  phases 
have  been  observed  in  the  chick  wing  and  other  structures 
(Wolpert,  1978;  Wilby  and  Ede,  1976). 

6.2.6  Summary  of  Experimentation  with  Growth 

The  processing  experiment  has  established  that  there  is 
no  apparent  distortion  in  distribution  of  cell  lineages 
(clones)  due  to  processing  order;  to  this  extent  the 
sequential  lumped  model  for  the  parallel  base  model  of 
growth  is  valid.  Further  validation  will  depend  on  further 
input-output  observation  in  the  context  of  intercellular 
processes.  Some  guidelines  for  time  step  management  have 
been  developed.  Options  for  growth  in  contiguous  structures 
have  been  described,  as  have  the  options  of  clocked  cell 
division,  division  in  response  to  the  threshold  of  a 
function,  and  the  Alberghina  division  decision  model. 
Possible  mechanisms  for  the  nonlocal  control  of  shape  and 
size  have  been  discussed,  and  an  example  given  of  shape 
mediated  by  direction  of  cell  division  under  the  influence 
of  i nterce 1 1 u 1 ar  communication. 

It  should  be  noted  that  these  simulation  results  cannot 
be  compared  directly  with  others,  since  they  are  unique  with 
respect  to  the  combination  of  individual  cellular  response 
( intracel lular  modelling)  and  changes  in  structures  of  cells 
( intercel lular  modelling).  The  cell  space  models  of  Wilby 
and  Ede  (1976)  and  Ransom  (1977)  described  in  Section  3.5, 
provide  for  individual  decisions  to  divide,  but  influences 
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on  such  divisions  cannot  be  locally  mediated.  Models  such 
as  that  of  Alberghina  (1977),  described  in  Section  3.5  and 
above,  describe  intracellular  (local)  influences  but  are  not 
applied  to  a  structure  of  cells;  their  domain  is  statistical 
measures  on  a  cell  population.  The  modelling  options 
presented  here  are  unique  in  providing  the  opportunity  to 
combine  both  levels  of  modelling. 

An  interesting  sideline  of  the  development  of  the 
growth  model  is  that  the  growth  of  cancerous  cells'  can  be 
studied.  A  cancerous  cell  is  one  in  which  division  is  no 
longer  properly  controlled  (Kim,  1975);  in  the  Alberghina 
model,  this  means  having  a  doubling  time  smaller  than  the 
"background"  cells.  Figure  6.10  illustrates  the  accelerated 
division  of  cell  lineage  "7"  in  contrast  to  normal  cells. 
Cell  lineage  7  has  a  doubling  time  of  1,  while  background 
cells  have  a  doubling  time  of  8  time  steps.  A  cancer 
researcher  might  experiment  with  sizes  of  tumours  as  related 
to  the  relative  doubling  times  of  normal  and  cancerous 
cells,  and  relate  his  results  to  an  intracellular  model  of 
di vi sion  control . 

6.3  Modelling  Interce 1 1 u 1 ar  Communication 

The  accumulating  evidence  for  positional  specificity  of 
cells  in  developmental  systems  necessitates  a  search  for 
processes  in  such  systems  capable  of  developing  local  values 
for  a  cell  which  define  it  relative  to  other  cells. 
Chemical  species  involved  in  the  competitive  activities  of 
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FIGURE  6.10  Lineage  7  divides  without  control 
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reaction  and  diffusion  have  been  considered  in  Section  3.4 
as  mechanisms  capable  of  producing  positional  specificity 
(pattern)  from  an  initially  uniform  stable  configuration. 

A  naive  attempt  to  model  the  transmission  of  signals 
between  cells  was  initially  made,  in  which  a  proportion  of 
the  current  concentration  of  a  cell  was  shared  between  its 
neighbours  at  each  time  step.  Whether  each  cell  possessed 
in  addition  a  reactive  capacity  to  produce  the  species  or 
not,  the  process  inevitably  produced  exponentially 
increasing  concentrations.  The  difference  equations 
representing  this  process  are  unstable  except  for 
prohibitively  small  time  steps.  They  may  be  regarded  as 
approximations  to  the  diffusion  equation 

u' (r,t)=L(u(r,t) ) , 

where  u'  is  the  partial  derivative  with  respect  to  time  and 
L(u)  is  the  Laplacian  operator  with  respect  to  the  spatial 
coordinates  r.  The  stability  criterion  for  explicit 
approximations  of  such  an  equation  is  that  the  time  step  be 
proportional  to  the  square  of  the  spatial  step,  which  makes 
explicit  approximation  impractical.  Instability  for  larger 
time  steps  yields  meaningless  results  as  numerical  errors 
increase  exponentially  with  time. 

Implicit  methods  for  solving  partial  differential 
equations  involve  adopting  a  difference  approximation  to  the 
differential  equation  in  which  the  unknown  quantity, 
u(r,t+1),  appears  explicitly  on  the  right  of  the  equation. 
For  each  point  in  a  spatial  difference  grid  such  an  equation 
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is  specified,  and  the  resulting  system  of  simultaneous 
equations  is  solved  at  each  time  step.  The  stability  of 
such  a  method  depends  on  the  eigenvalues  of  the  matrix  of 
coefficients  of  u(r,t+1);  if  no  eigenvalue  exceeds  1  the 
process  is  stable  (Fox,  1962). 

Reaction-diffusion  systems  of  interest  are  likely  to  be 
nonlinear  (Grierer  and  Meinhardt,  1972)  and  coupled  (Gmitro 
and  Scriven,  1966,  Kauffman,  1978).  Nonlinearity  may  be 
treated  in  finite  difference  approximations  by  iterative 
methods  of  solution  of  the  resulting  nonlinear  difference 
equations;  however,  the  analysis  of  departures  from  stable 
equilibria  make  the  use  of  quas i 1 i near i zed  forms  of  the 
equations  appropriate  (Gmitro  and  Scriven,  1966). 

Coupling  of  equations  may  be  dealt  with  either  by 
representing  the  coupling  terms  explicitly  in  the  linear 
difference  equations  or  by  the  use  of  predi ctor -corrector 
methods.  The  former  approach  necessitates  m2n2  storage, 
where  m  is  the  number  of  coupled  species  and  n  the  number  of 
points,  and  mnt  time  units,  where  t  is  the  average  time  to 
solve  for  one  unknown  if  an  elimination  method  of  solution 
is  used.  The  predictor-corrector  approach  uses  much  less 
space,  mn2,  but  the  time  for  solution  is  kmnt  where  k  is  the 
number  of  uses  of  predictor  and  corrector  until  convergence 
is  obtained  at  each  time  step.  Since  the  stability 
criterion  is  such  that  iterative  methods  of  solution  of  the 
system  of  equations  are  effective,  time  is  less  at  a  premium 
than  space  and  the  predi ctor -cor rector  approach  is  selected. 
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The  Crank-N i cho 1  son  approximation  (Fox,  1962)  was 
judged  appropriate  for  the  solution  of  the  linearized  system 
A: 

ul'  (r,t)=K( 1 , 1 ) u 1 +K ( 1 , 2 ) u2+ .  . . +K ( 1 ,n)uN  +  D( 1 ) L ( u 1 ) 
u2' (r,t)=K(2, 1 )u1+K(2,2)u2+. . . +K ( 2 , n ) uN+d ( 2 ) L ( u2 ) 


uN'  (r,t)=K(N,1)u1  +  ..  . +K ( n , n ) uN+d ( N ) L ( uN ) 

In  the  approximation,  the  quantities  on  the  right  hand 
side  are  replaced  by  the  average  of  their  values  at  t+1  and 
t,  the  Laplacian  is  replaced  by  central  difference  formulae 
for  the  second  derivatives,  and  the  left  hand  side  is 
replaced  by  a  single  forward  difference  approximation.  The 
method  has  local  truncation  error  0(k2+h2)  in  regions  free 
of  singularities  and  therefore  combines  accuracy  and 
s tabi 1 i ty . 

Boundary  values  are  incorporated  into  the  linear  system 
of  equations,  care  being  taken  that  the  eigenvalues  of  the 
resulting  matrix  not  exceed  zero.  Boundary  values  in 
problems  of  interest  are  either  normal  flux  equal  to  zero  at 
the  boundary  (no  communication  with  outside  regions)  or 
specified  values  of  concentrations  on  the  boundary. 
Together  with  initial  values  for  u(r,t)  at  the  initial  time, 
they  ensure  a  unique  analytic  solution  to  the  differential 
system,  and  therefore,  under  conditions  of  convergence  and 
stability,  to  the  difference  system. 
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6.3.1  Formal  Specification  of  the  Communication  Model 

A  communication  model  is  a  structure 

S= <Q , T , i , W , d> ,  where 

Q  is  the  single  state  variable,  consisting  of  a  linked  list 
of  cells  specified  by  neighbourhood  links  and  pointers  to 
individual  parameters; 

T  is  the  continuous  time  base; 

I  is  the  input  value  set,  incorporating  communication 
through  the  boundaries  of  the  structure  with  the  environment 
(boundary  values); 

W  is  the  input  segment  set,  i ncorpor at i ng  the  possible 

dependence  of  the  boundary  values  on  time; 
d  is  the  state  transition  function,  d:Q  x  I  -->  Q 

The  state  transition  function  d  generates  the  new  state 
q  by  changing  one  or  more  of  the  parameters  pointed  to  by 
each  cell.  Its  domain  is  the  entire  structure,  since  change 
in  each  cell  depends  on  the  state  of  neighbouring  cells.  In 
the  present  instance,  d  is  given  in  differential  form  in  a 
reaction-diffusion  system  for  the  concentration  of  one  or 
more  chemical  species. 

The  continuous  time  system  S  can  be  simulated  exactly 
by  a  discrete  time  system  SD  provided  that  its  state 
trajectories  are  to  be  reproduced  only  at  sample  intervals 
of  constant  length  h  (Ziegler,  1976).  However  the 
trajectories  can  be  generated  only  if  the  global  transition 
function  (that  is,  the  analytic  solution  of  the  differential 
system)  is  known.  SD  is  therefore  unobtainable  as  a 
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practical  basis  for  iterative  generation  of  state 
tra jector i es ,  but  may  be  used  as  a  standard  for  judging 
approximate  methods. 

The  system  SD  is  given  by 
SD= <Q , TD , I , Wh , dh> ,  where 
Q,  I  are  as  in  S; 

TD  is  a  discrete  time  base; 

Wh  is  the  input  segment  set  W  restricted  to  segments  of 
length  h; 

dh  is  the  state  transition  function  operating  at  the 
endpoints  of  input  segments  (at  discrete  time  steps  nh ) . 

Correspondi ng  to  the  system  SD  there  is  an  iterative 
specification  (Ziegler,  1976),  and  correspondi ng  to  the 
iterative  specification  a  sequential  machine  Md  which 
simulates  it  exactly  (generates  the  state  trajectories  given 
the  global  transition  function).  Let  MD(h)  be  given  by 
<Wh,Q,dh>.  Now  an  approximation  method,  which  represents 
the  differential  system  as  a  system  of  difference  equations, 
may  be  referred  to  as 

MA ( h ) =< I , Q , da ( h ) > ,  where  the  transition  function 
da(h):Q  x  I  -->  Q 

is  given  by  the  solution  of  the  difference  system  for 
q(  t+h)  . 

It  is  necessary  to  establish  an  "approximate" 
specification  morphism  between  MA(h)  and  MD(h)  in  order  that 
the  lumped  numerical  model  MA(h)  may  be  regarded  as  a  valid 
simulator  of  the  lumped  model  MD(h),  and  therefore  of  the 
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base  model  S  in  a  particular  experimental  frame.  A 
specification  morphism  between  the  two  sequential  machines 
( preservat i on  of  transition  function)  is  established  by 
proving  that  the  trajectories  generated  by  the  approximate 
machine  MA(h)  come  arbitrarily  close  to  those  generated  by 
MD(h)  on  identical  input  segments. 

Conditions  under  which  this  situation,  Known  as 
"convergence"  in  numerical  mathematics,  will  hold  true 
depend 'on  "accuracy"  and  "stability".  The  accuracy  of  a 
numerical  method  depends  on  its  local  truncation  error;  that 
i  s  ,  on 

| dh (g , w) -da ( h ) (g , w( 0 ) |=M 

In  the  case  of  the  Crank-Ni chol son  method,  the 

truncation  error  is  0(h2+K2).  The  stability  of  a  numerical 
method  depends  on  the  propagation  of  this  error,  and  the 

errors  incurred  by  the  use  of  finite  mathematics,  by 

repeated  application  of  the  method.  In  the  Crank-Ni col  son 
approximation  the  eigenvalues  of  the  matrix  of  coefficients 
in  the  difference  equations  are  less  than  one,  which  implies 
that  repeated  solution  of  the  system  does  not  allow 

unbounded  error  amplification  (Fox,  1962).  In  terms  of 
MD(h)  and  MA(h)  this  means  the  accumulated  error  after  n 
steps 

| dh ( g , w 1 , w2 , . . . , wn ) -da ( h ) ( g , wl ( 0 ) , w2 ( 0 ) , . . . , wN ( 0 ) ) | 
is  bounded  by  Mp(h)  where  p  is  a  polynomial  expression  in  h. 
Convergence  is  then  assured  as  Mp(h)  -->  0  as  h  -->  0,  and 
the  error  is  arbitrarily  small. 


* 

B9 


■ 


.. 1 . 


160 


Thus,  because  the  characteristics  of  the  Crank- 
Nicholson  approximation  are  Known,  it  may  be  asserted  that 
an  approximate  specification  morphism  holds  between  MA(h) 
and  MD ( h )  (between  the  numerical,  iterative  solutions  of  the 
difference  equations  and  the  analytic  solutions  of  the 
differential  equations)  in  any  experimental  frame. 

The  real  significance  of  this  fact  is  that  experiments 
run  by  the  simulator  generating  trajectories  may  be 
considered  direct  tests  of  the  base  model  (reaction- 

diffusion)  in  comparison  with  trajectories  of  the  real 
system  in  a  particular  experimental  frame;  that  is,  the 
effects  of  replacing  the  differential  system  by  a  difference 
(spatially  and  temporally  discrete)  system  may  be  ignored. 
Experiments  with  reaction-diffusion  as  the  mechanism  of 
communication  between  cells  may  therefore  concentrate  on 
observing  generated  trajectories  (  a  sequence  of  "patterns") 
and  comparing  them  with  events  in  real  systems,  thus 
investigating  the  equivalence  of  their  input-output 
relations  (postulate  10,  Ziegler  (1976)). 

It  may  be  noted  that  the  sequential  machine  MA(h)  is 
identified  directly  with  the  simulator  or  computer  program. 
MA(h)  operates  in  parallel  on  the  cell  structure,  updating 
all  cells  simultaneously  for  each  time  step.  The  implicit 
equations  are  solved  as  a  simultaneous  system  in  the 
computer  program.  Parallel  processing  is  thus  preserved  and 
the  identification  is  justified. 

The  communication  model  is  specified  only  at  the  level 
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of  the  cell  structure,  in  contrast  to  the  growth  model. 

6.3.2  Experiments  with  Numerical  Solutions  of  Reaction- 
Diffusion  Systems 

The  concept  of  positional  specificity  in  pre¬ 
patterning  cells  for  differentiation  is  a  hypothesis 
inferred  from  many  events  in  development,  as  discussed  in 
Chapter  Three.  It  has  proved  impractical  to  obtain  direct 
evidence,  in  the  form  of  local  readings  on  concentration 
values,  of  such  fields  due  to  the  multitude  of  processes 
occurring  in  cells  and  the  difficulty  of  performing  in  vivo 
experiments.  The  "real  system",  therefore,  cannot  be  given 
by  concentration  data  from  a  functioning  morphogenetic 
field.  For  the  purposes  of  demonstrating  the  efficacy  of 
reaction-diffusion  as  a  mechanism  for  generating  positional 
specificity,  the  real  system  will  be  considered  to  be  any 
group  of  cells  displaying  such  specificity,  that  is, 
nonuni formi ty . 

Prepatterning  hypothesis:  Reaction-diffusion  is  a  mechanism 
capable  of  producing  stable  nonuniform  configurations  from 
uniform  configurations. 


for  experimentation  with 


In  order  to  provide  a  background 
the  prepat terni ng  hypothesis 


a  brief  review  of  the 
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instability  analysis  for  R-D  systems  is  justified.  Reaction 
and  diffusion  are  "competing"  processes  in  that  reaction  is 
disruptive  (produces  local  change)  whereas  diffusion  is 
dissipative  (distributes  i r regu 1 ar i t i es ) .  This  interaction 
may  be  expressed  in  nonlinear  (coupled)  differential 
equations  by  means  of  the  relative  magnitude  of  the  reaction 
rates  and  diffusion  rates.  The  solution  of  such  equations 
may  be  expressed  as  the  sum  of  a  uniform,  steady  state 
solution  and  an  excursion  from  that  state.  The  differential 
system  may  then  be  adjusted  to  describing  the  excursion  from 
steady  state  by  means  of  the  process  of  quas i 1 i near i zat i on 
about  the  steady  state  (Gmitro  and  Scriven,  1966).  The 
spatial  dependence  of  the  solution  for  the  excursion  in  a 
domain  may  be  expressed  as  the  sum  of  the  Laplacian  shape 
functions  relevant  to  that  domain,  multiplied  by  an 
exponential  growth  factor  related  to  the  relative  magnitudes 
of  the  reaction  and  diffusion  parameters  (Gmitro  and 
Scriven,  1966,  Kauffman  et  al,  1978).  In  general  this  sum 
is  zero  and  the  steady  state  is  maintained;  however,  at 
critical  sizes  of  the  domain,  one  particular  shape  function 
may  fit  exactly  into  the  domain;  it  is  then  selected  for 
amplification  and  causes  a  departure  from  uniform  steady 
state . 

Gmitro  and  Scriven  (1966)  discuss  qualitatively 
criteria  under  which  such  events  occur  in  systems  of  one, 
two  or  three  participating  species.  For  N  participating 
species  such  behavior  can  be  ensured  by  only  (N*N+3*N)/2 
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constraints  on  2N*N+1  parameters.  Thus  there  is  a 
considerable  degree  of  freedom  in  selecting  values  for  the 
reaction  and  diffusion  parameters.  In  terms  of  the  real 
system  this  implies  versatility,  but  in  terms  of 
experimentation  it  implies  the  impr act i ca 1 i ty  of  exhaustive 
investigation  of  behaviors  of  reaction-diffusion  systems. 
However,  for  the  purposes  of  testing  the  prepat terni ng 
hypothesis,  it  is  sufficient  to  establish  such  behavior  in 
only  one  system.  The  system  chosen  for  investigation  is 
that  of  Kauffman  et  al.  (1978)  given  by 
u 1 '  =K ( 1 , 1 )u1+K( 1 , 2 ) u2  +  D ( 1 ) L ( u 1 ) 
u2'  =K ( 2 , 1 )u1+k(2,2)u2+D(2)L(u2) 
where  K( 1 , 1 )=-7.8,  K(1,2)=13.4,  K ( 2 ,  1  )  = -  1 ,  K(2,2)=2.04 
D  (  1  )  =  1 6  ,  D  (  2  )  =  1  . 

Analysis  by  the  methods  of  Gmitro  and  Scriven  (1966) 
predicts  that  the  critical  size  range  is  the  interval 
(6. 1,9.1)  (Kauffman  et  al.  1978). 

Prepat terni ng  experiment:  To  establish  the  prepat terni ng 
hypothesis  by  observing  departures  from  uniform  steady  state 
in  a  reaction  diffusion  system,  namely  the  Kauffman  system 
given  above. 

Experimental  frame:  The  Kauffman  system  is  to  be  applied  in 
one-  and  two-dimensional  rectangular  domains.  The  output 
trajectories  for  various  sizes  are  to  be  compared  to  the 
expected  shape  functions.  The  boundary  conditions  are  that 
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the  normal  flux  be  zero,  corresponding  to  an  autonomous 
system  (no  communication  through  the  boundary);  the  shape 
functions  expected  for  rectangular  domains  under  these 
boundary  conditions  are  integral  numbers  of  half-cosine- 
waves  (Kauffman  et  al.,  1978). 

Procedure:  The  Crank-Nicol son  method  is  used  in  generating 
trajectories  in  one-  and  two-dimensional  domains.  In  the 
case  of  one  dimension,  successive  critical  sizes  are 
investigated  by  manipulating  h,  here  representing  the 
spatial  step  size  or  distance  between  cell  locations,  and  n, 
the  number  of  cells,  bearing  in  mind  that  h  must  be  less 
than  1  to  ensure  accuracy.  In  the  case  of  two  dimensions, 
since  the  Crank-Nicolson  method,  with  its  stability 
character i st i cs ,  is  represented  in  a  rectangular  coordinate 
system,  one  of  two  actions  is  necessary.  Either  the  Crank- 
Nicholson  method  must  be  rephrased  in  terms  of  a  hexagonal 
array,  or  the  hexagonal  array  must  be  converted  to  a 
rectangular  one.  The  second  alternative  is  chosen,  firstly, 
because  the  Laplacian  operator  in  hexagonal  coordinates 
exhibits  mixed  partial  derivatives  which  prohibit  the  use  of 
alternate-direction  acceleration  techniques  (Fox,  1962) 
which  may  prove  useful  in  some  circumstances,  and  secondly, 
because  there  is  no  guarantee  of  stability  in  the  resulting 
coefficient  matrix.  Conversion  of  the  hexagonal ly  arranged 
cells  into  a  rectangular  array  for  the  purposes  of  reaction- 
diffusion  in  some  of  their  parameters  necessitates  the 
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initial  interpolation  of  extra  values  for  these  parameters 
at  intermediate  points.  Aitken  interpolation  in  two 
directions  is  used  to  achieve  this,  care  being  taken  that 
enough  points  are  used  to  maintain  the  accuracy  of  the 
Crank-Ni chol son  method.  When  a  trajectory  is  completed,  the 
extra  values  are  discarded  and  the  concentration  values 
returned  to  the  parameter  records  of  the  correspondi ng 
cells. 

Results:  The  expected  patterns  were  observed  in  the  expected 
size  ranges,  and  are  presented  in  Figures  6.11-6.14. 

Discussion:  The  trajectories  generated  by  the  reaction 
diffusion  equation  display  a  dependence  on  various 
parameters,  both  of  the  numerical  approximation  and  of  the 
differential  system  itself. 

Refinements  of  the  spatial  grid  and  of  the  time  step 
produce  faster  convergence  of  the  predictor -corrector  pair 
used  to  represent  the  coupling  of  the  two-species  equation, 
but  increase  the  order  of  the  system  to  be  solved  and  number 
of  time  steps  to  be  taken;  an  optimal  balance  has  to  be 
struck,  minimizing  the  number  of  arithmetic  operations.  The 
unconditional  stability  of  the  method  is  demonstrated  by  the 
freedom  with  which  such  manipulation  can  be  performed  with 
regard  only  to  the  accuracy  constraint  that  the  steps  be 
less  than  1 . 

The  speed  with  which  the  patterns  reach  within  0.1  of 
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FIGURE  6.11a  First  cosine  pattern  arises  in  one 
dimensional,  two  species  system 
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SPECIES  1  PLOTTED  FOE  1  OOT  OF  1 

CONCENTRATION  SCALED  BY  1.000000 

SECTION  AL0N3  X  AXIS 


FIGURE  6.11b  First  cosine  pattern  established,  species  one 
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SPECIES  2  PLOTTED  FOR  1  OUT  OF  1  CELLS 
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FIGURE  6.11c  First  cosine  pattern,  species  two 
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SPECIES  1  PLOTTED  FOR  1  OUT  OF  1  CELLS 

CONCENTRATION  SCALED  BY  1,000000 

SECTION  ALONG  X  AXIS 


FIGURE  6.12a  Second  cosine  pattern,  species  one 
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FIGURE  6.12b  Second  cosine  pattern,  species  two 
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FIGURE  6.13  Third  cosine  pattern,  species  one 
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SPECIES  1  PLOTTED  FOP  1  ODT  OF 
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FIGURE  6.14  Fourth  cosine  pattern,  species  one 
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the  pure  cosine  functions  varies  with  the  initial 
disturbance  (noise)  induced  locally  in  the  uniform  steady 
state  to  provide  the  impetus  for  departure  from  uniformity. 
It  may  be  noted  that  when  noise  is  not  modelled  directly, 
the  accumulating  numerical  error  is  sufficient  eventually  to 
initiate  the  formation  of  the  non-uniform  pattern.  The 
system's  effect  on  numerical  error  is  interesting  in 
general.  In  uniform  ( noncr i t i ca 1 )  states,  the  Laplacian 
term  "distributes"  the  error  in  the  system,  thus  dissipating 
it;  at  critical  sizes,  numerical  error  acts  constructively 
in  establishing  the  pattern,  and  then  in  the  new  uniform 
configuration  is  again  dissipated.  A  speculation  may  be 
advanced  as  to  whether  the  discrete  representat i on  of 
reaction-diffusion  perhaps  has  a  structural  validity  for  the 
real  system  (that  is,  that  cells  should  be  regarded  as 
discrete  entities  and  not  as  points  in  a  continuous 
di f f erent i a  1  field). 

A  disturbance  in  only  one  cell  of  a  field  is  sufficient 
to  induce  the  pattern.  Gmitro  and  Scriven  (1966)  suggest 
that  noise  in  the  real  system  is  caused  by  "thermal 
fluctuations"  in  the  concentrations  due  to  constituent 
molecule  Kinetics.  In  a  developing  system,  however,  another 
source  of  variation  is  provided  by  dividing  cells.  During 
mitosis,  copies  are  made  of  some  substances  while  others  are 
distributed,  perhaps  not  evenly,  between  daughter  cells.  A 
local  disturbance  in  concentration  results.  It  may  well  be 
that  cell  division,  taking  a  morphogenetic  field  from  one 
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pattern  to  the  next  as  the  size  increases,  also  controls  the 
speed  with  which  patterns  arise  and  decay  by  providing  noise 
of  a  certain  magnitude.  If  so,  it  is  another  instance  of 
the  essential  interdependence  of  growth  and  determination  in 
developmental  systems. 

The  system  explored  in  the  prepatterning  experiment  may 
be  thought  to  be  highly  artificial;  in  reality  there  is  such 
a  high  degree  of  freedom  in  the  system  that  it  may  be  used 
to  represent  many  actual  configurations.  Firstly,  critical 
lengths  may  be  achieved  in  fields  of  any  desired  size  by 
varying  the  step  size  (distance  between  cells)  in  inverse 
proportion  to  the  number  of  cells.  Secondly,  the  K  and  D 
parameters  are  chosen  to  satisfy  certain  conditions  for  the 
desired  stability  behavior  (Gmitro  and  Scriven,  1966),  such 
conditions  in  this  case  imposing  5  constraints  on  6 
quantities.  The  remaining  degree  of  freedom  permits  the 
system  parameters  to  be  scaled  downwards  to  values  realistic 
for  cell  reactions  and  diffusions.  The  experimenter  may 
take  one  of  two  approaches;  he  may  postulate  values  for  the 
parameters  in  the  study  of  a  particular  real  morphogenic 
field,  and  see  whether  the  resultant  patterns  are 
established  after  a  realistic  time  period,  or  he  may  use 
knowledge  about  the  actual  time  taken  in  establishing  a 
field  in  order  to  "calibrate"  the  system  parameters. 

While  the  Kauffman  system  cannot,  obviously,  represent 
all  stability  behaviors  of  interest  in  pattern  formation 
models,  it  does,  because  of  the  discussed  degrees  of 
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freedom,  have  the  capability  to  represent  (at  least 
qualitatively)  events  occurring  in  many  developing  systems. 
This  property  of  the  system  is  used  to  advantage  in  Section 
6.5. 


Conclusions:  Experimentation  has  established  the 
prepat terni ng  hypothesis,  that  is,  reaction-diffusion  is  a 
mechanism  capable  of  producing  nonuniform  configurations 
from  un i formi ty . 

The  use  of  numerical  techniques  becomes  costly  with 
many  points,  particularly  in  two-dimensional  systems  where 
there  i s  translation  to  a  rectangular  grid.  The  question 
arises  as  to  whether,  in  domains  whose  shape  functions  are 
Known,  iterative  generation  of  the  trajectories  in  the 
growing  domain  is  necessary.  Under  certain  conditions  it 
is.  If  the  appearance  and  decay  of  the  patterns  are  to  be 
studied  for  themselves,  or  if  other  processes  being  modelled 
depend  dynamically  on  local  concentrations,  then  iterative 
generation  is  essential. 


6.3.3  Modelling  Using  Analytical  Solutions 

In  some  modelling  situations  the  communication  model 
may  be  independent  of  any  other  model  and  the  shape 
functions  for  the  system  may  be  available.  The  reaction- 
diffusion  continuous  system  S  may  then  be  modelled  by  the 
machine  MD(h)  described  in  Section  6.3.1,  the  exact  solution 
machine  giving  trajectories  at  discrete  values  of  time.  The 
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time  base  of  MD(h)  has  variable  step  length,  the  length  of 
any  one  step  being  given  in  the  continuous  system  by  the 
time  between  establishment  of  one  configuration  and  the 
next.  An  "event"  in  MD ( h )  is  initiated  when  critical  size 
is  reached,  that  is,  when  a  dynamically  mediated  condition 
on  the  structure  of  the  cell  system  is  satisfied. 

The  option  of  using  analytic  solutions  for  reaction- 
diffusion  systems  in  one-dimensional  and  two-dimensional 
rectangular  domains  is  provided.  The  relevant  shape 
functions  are  sines  or  cosines,  as  is  demonstrated  in 
Figures  6.15  and  6.16.  In  Figure  6.15a  and  6.15b,  the 
solutions  for  one  of  the  species  in  orthogonal  directions 
are  given,  and  in  Figure  6.15c  the  result  of  multiplying 
these  components  is  shown  by  indicating  the  sign  at 
individual  cells  in  the  two-dimensional  structure.  Since 
these  patterns  of  signs  at  individual  cell  indicate  the 
location  of  lines  where  the  solutions  are  zero,  they  are 
referred  to  as  "nodal  patterns".  Similarly  Figures  6.16a-c 
show  the  solution  for  another  critical  size. 

The  shape  functions  for  circular  and  elliptic  domains 
are  Bessel  and  Matthieu  functions  respectively  (Kauffman  et 
a  1 . y  1978).  The  nodal  patterns  which  may  arise  given 
certain  critical  sizes  of  the  radii,  (or  the  major  and  minor 
axes),  are  Known  qualitatively  (diameters  and  concentric 
circles,  and  axes,  hyperbolae  and  confocal  ellipses),  but 
finite  expressions  for  evaluating  them,  thus  locating  and 
timing  the  rise  of  nodal  lines,  are  difficult  to  obtain. 
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FIGURE  6.15a  Analytic  solution,  x  direction  (pattern  one) 
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FIGURE  6.15b  Analytic  solution,  z  direction  ( noncr i t i ca 1 ) 
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FIGURE  6.15c  Nodal  pattern 
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FIGURE  6.16a  Analytic  solution,  x  direction  (pattern  one) 
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FIGURE  6.16b  Analytic  solution,  z  direction  (pattern  one) 
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Figure  6.16c  Nodal  pattern 
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Sometimes  a  nodal  configuration  may  be  inferred  by 
extrapol at  ion  from  rectangular  configurations  or  by  symmetry 
considerations,  as  will  be  seen  in  Section  6.5.  In  a  model 
where  the  positional  specificity  of  cells  is  specified 
completely  by  these  nodal  lines,  analytic  expressions  for 
the  shape  functions  themselves  are  not  necessary,  merely 
knowledge  of  their  sign. 


6.3.4  Summary  of  Experiments  with  Communication 

The  numerical  solution  of  reaction-diffusion  systems 
has  been  investigated  for  a  system  of  wide  applicability. 
The  method  generates  trajectories  approaching  those  of  the 
base  model  in  all  chosen  experimental  frames,  reflecting  the 
approximate  specification  morphism  between  the  discrete 
approximation  and  the  analytical  solution.  The  validity  of 
the  reaction-diffusion  model  for  producing  prepatterns  in 
morphgenetic  fields  must  be  tested  in  the  context  of 
observable  developmental  events  (for  example,  cell 
differentiation)  in  models  which  address  the  problem  of  cell 
determination  in  response  to  pre-patterns. 

The  options  of  the  use  of  analytic  solutions  or  node 
location  are  provided  for  certain  domains  and  situations. 


6.4  Modelling  Growth,  Communication  and 
Development  is  a  process  in 
communication  between  cells  operate 


Determi nation 
which  growth 
s imu 1 taneous ly , 


and 

and 


critical  conditions  (presumably  local)  at  various  stages 
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the  process  initiate  changes  in  the  organism  resulting  in 
differentiated  cell  types.  In  order  to  model  development, 
the  models  of  growth  and  communication  must  be  merged,  and  a 
model  developed  for  recognizing  critical  conditions. 

The  base  models  for  growth  and  communication,  G  and  S, 
have  continuous  time  bases  and  the  same  state  set,  the  set 
of  cell  structures.  It  is  therefore  a  simple  matter  to 
merge  them  into  one  system.  However  for  purposes  of 
simulation  their  discrete  lumped  models  in  particular 
experimental  frames  must  be  considered.  In  the  same 
experimental  frame,  the  time  bases  for  the  two  lumped  models 
must  be  constructed  so  as  to  ensure  that  possible  influences 
of  one  process  on  the  other  are  taken  into  account. 

In  Section  6.2  the  possible  influence  of  local 
conditions  on  a  cell's  capacity  to  divide  was  discussed. 
The  discrete  representat ion  of  the  growth  time  base  is  valid 
if  local  conditions  affecting  growth  do  not  change  during  a 
growth  time  step.  That  is  to  say,  if  G  depends  on  S  then 
timestep(G)  is  less  than  or  equal  to  timestep(S),  and 
timestep(S)  must  be  an  integral  multiple  of  timestep(G). 

The  effect  that  growth  has  on  the  communication  process 
may  be  negligible  or  profound.  If  one  or  two  extra  cells  in 
a  structure  do  not  take  the  system  into  any  new  critical 
size  range,  the  spatial  patterns  will  not  change;  however, 
if  the  expansion  does,  a  completely  new  configuration  will 
arise  as  concentration  levels  adjust  by  means  of  chemical 
flow  through  the  system.  Since  the  communication  process 
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must  be  sensitive  to  the  possibility  of  entering  a  new 
critical  size  at  ail  times,  S  depends  on  G.  Timestep(S) 
should  be  less  than  or  equal  to  timestep(G)  and  timestep  (G) 
must  be  an  integral  multiple  of  timestep  (S). 

Therefore  if  both  processes  are  dependent  on  each 
other,  they  must  have  the  same  discrete  time  step.  In  cases 
where  growth  is  independent  of  the  communication  process, 
communication  time  steps  can  be  smaller  than  growth  time 
steps. 

Determination,  the  genetic  commitment  of  a  cell  to  a 
future  path  leading  to  differentiation,  is  a  response  to 
attainment  of  critical  values  in  a  cell's  local  conditions. 
Change  in  such  conditions  is  generated  by  the  communication 
model.  A  model  for  determination,  therefore,  being 
sensitive  to  any  changes  as  they  arise,  should  have  the  same 
time  step  as  the  communication  model. 

6.4.1  Formal  Specification  of  the  Determination  Model 
A  determination  model  is  a  structure 
D=<Q,T,e>  ,  where 

Q  is  the  single  element  state  set,  consisting  of  a  linked 
list  of  cells  specified  by  neighbour  links  and  pointers  to 
individual  parameters; 

T  is  the  continuous  time  base; 
e  is  the  state  transition  function, 
e : Q  -->  Q. 

The  system  D  has  state  transition  function  e  which  (a) 
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leaves  the  structure  unchanged  if  no  (local)  critical 
parameters  are  found,  and  (b)  changes  the  character  of  cells 
in  the  structure  for  which  critical  parameters  are  found  (by 
manipulating  the  cell's  parameter  record,  for  example  giving 
it  a  new  "genetic  activation  code"). 

It  may  be  remarked  that  determination,  once  initiated, 
is  a  purely  local  event,  occurring  in  individual  cells  as  an 
isolated  response,  neither  changing  the  structural 

arrangement  of  the  group  of  cells  nor  changing  i nterce 1 1 u 1 ar 
communication.  Providing  the  time  base  is  such  that  local 
conditions  do  not  change  during  the  process  of 
determination,  it  is  immaterial  whether  the  cells  are 
treated  sequentially  or  in  parallel  by  the  transition 
function  e.  D  may  therefore  be  modelled  by 
DS  =  <Q,TD,ed>,  where 
TD  is  a  discrete  time  base; 
ed  is  the  state  transition  function; 

ed : Q  -->  Q,  where  ed  is  given  by  the  following  algorithm: 

1.  Set  pointer  P  to  the  first  cell  on  the  processing  list. 

2.  Inspect  relevant  parameters  of  the  cell  for  critical 
va 1 ues . 

3.  If  a  critical  value  is  found,  effect  the  resulting 
changes  in  the  cell's  character. 

4.  Set  P  to  the  next  cell  on  the  processing  list. 

5.  If  P  is  not  null,  go  to  2,  otherwise  stop. 

The  particular  changes  in  a  cell's  character  brought 
about  by  determination  will  depend  upon  the  experimental 
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frame  within  which  a  particular  experiment  with  DS  is  being 
undertaken.  It  may  take  the  form  of  adding  digits  to  the 
developmental  genetic  code  (see  Section  4.1),  or  noticing  a 
significant  accumulation  of  such  digits  and  initiating 
actual  differentiation;  the  latter  situation  may  result  in 
division  of  the  structure  of  cells  into  substructures ,  which 
possibly  do  not  communicate  similarly.  The  reaction- 
diffusion  system  in  progress  may  be  eliminated  in  the 
differentiated  cells.  Thus  a  still  higher  degree  of  control 
is  implied  in  a  model  of  differentiation.  In  Section  6.5 
such  systems  will  be  discussed. 

6.4.2  Growth  and  Communication 

An  experiment  has  been  conducted  with  one-dimensional 
structures  in  which  cell  division  by  the  Alberghina  model 
occurs  simultaneously  with  communication  using  the  Kauffman 
system.  Two  situations  are  considered:  (a)  the  patterns 
arising  are  allowed  to  stabilize  before  further  cell 
division  is  permitted  (timestep  (S)  much  less  than  timestep 
( G ) ) ,  and  (b)  cell  division  is  allowed  to  occur  during  the 
establishment  (and  decay)  of  patterns  (timestep  (S)  less 
than  timestep  ( G ) ) .  Representative  results  may  be  seen  in 
Figure  6.17,  where  a  one-dimensional  structure  is  growing 
and  undergoing  various  non-uniform  configurations.  It  is 
found  that  the  effect  of  dividing  cells  on  the  dynamics  of 
pattern  formation  is  as  expected  by  Section  6.3.2.  Cell 
division  (a)  helps  a  new  pattern  to  rise  by  providing  noise" 
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FIGURE  6.17a  Growth  through  cosine  patterns,  initial 
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FIGURE  6.17b  With  10  cells 
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FIGURE  6. 17c  With  14  cel  Is 
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Figure  6 .  1 7d  With  20  cells 
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FIGURE  6. 17e  With  40  cel  Is 
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(simulated  by  sharing  the  concentration  of  the  relevant 
substance  unevenly),  (b)  does  not  affect  the  time  taken  in 
establishing  the  pattern  (the  local  "flattening"  is 
compensated  for  by  the  additional  noise  supplied  to  the 
process)  unless  the  critical  size  for  the  system  is 
exceeded,  and  (c)  hastens  the  decay  of  patterns  by  local 
"flattening" . 

It  should  be  remarked  that  the  rate  of  growth  has  a 
direct  effect  on  the  patterns  appearing  in  a  system.  If  the 
system  grows  fast  it  may  miss  some  of  the  patterns  (for 
example,  by  going  from  10  to  20  cells  in  Figure  6.17,  thus 
missing  an  intermediate  pattern).  In  such  cases  patterns 
begin  to  arise  but  the  system  exceeds  the  critical  size 
before  the  pattern  stabilizes.  A  comparison  may  be  made 
between  the  capabilities  exhibited  in  Figure  6.17,  and  the 
capabilities  of  systems  like  CELIA  (Herman  et  al.,  1974) 
which  also  investigates  one-dimensional  structures  (see 
Section  3.4.1).  Firstly,  the  state  set  for  this  author's 
model  may  be  extensive  and  of  mixed  type,  while  the  state 
set  of  CELIA  has  a  single  component,  a  discrete  genetic 
marker.  Secondly,  in  CELIA  signal  transmission  is  of  a 
primitive  nature.  Thirdly,  due  to  this  primitive 
signalling,  gradient  patterns  are  achieved  only  by  imposing 
special  states  at  the  end  points  of  arrays,  thus  violating 
strict  initial  uniformity.  In  summary,  CELIA  presents  a 
much  less  "realistic"  picture  of  developing  one-dimensional 
systems . 
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6.5  Modelling  Developmental  Systems 

A  developmental  system  is  one  in  which  the  various 
intercellular  processes  interact  in  a  manner  resulting  in 
organized  changes  in  cells  and  cell  structures.  In  the 
context  of  this  project,  these  processes  have  been 
classified  under  the  headings  growth,  communication  and 
determination.  Development  of  the  system  appears  to  occur 
by  means  of  coordination  and  control  of  these  processes. 

It  would  be  easy  to  postulate  a  "global  control"  which 
performs  the  task  of  coordination,  making  sure  that  the 

growth  rate  takes  the  system  through  the  right  patterns, 

% 

scheduling  determination,  and  so  on.  In  other  words,  it 
would  be  easy  to  specify  a  master  system  invoking  growth, 
communication  and  determination  when  required.  However  to 
do  so  would  be  to  act  in  contradiction  of  the  basic 
assumption  that  development  is  mediated  at  the  local  level, 
though  its  effects  are  widespread. 

How  is  it,  then,  that  such  processes  act 
"cooperatively"  to  produce  a  viable  system?  The  answer  must 
lie  in  the  process  of  natural  selection:  many  nonviable 
interactions  must  have  been  discarded  in  selecting  viable 
systems  for  survival.  It  is  possible  to  imagine  the  genesis 
of  a  simple  organism,  for  example.  A  capacity  for  a  certain 
reaction-diffusion  system  is  inherited  or  provoked  by 
mutation  in  a  cell.  It  divides  and  happens  to  reach  a 
critical  shape  for  that  R-D  system.  The  resulting 
nonuniform  pattern  triggers  changes  in  some  of  the  cells. 
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This  configuration  is,  in  some  sense,  viable  in  its 
environment  and  therefore  is  selected  for  and  continues  to 
develop.  That  reaction-diffusion  system,  after  surviving 
many  such  viability  tests,  becomes  one  of  the  developmental 
tools  of  the  strain's  evolutionary  descendents. 

In  combining  all  of  the  processes  of  development  in 
order  to  model  a  particular  real  system,  then,  what  appears 
to  be  a  controlled,  cooperative  interaction  must  be  thought 
of,  rather,  in  terms  of  fortuitous  selection  of  system 
parameters.  The  growth  and  communication  time  bases  must  be 
selected  with  reference  to  the  particular  experimental 
frame,  and  i nterpretat ion  of  determination  decided  for  that 
context.  Thus  there  will  be  no  formal  specification  of  a 
developmental  model,  as  such  a  structure  bears  no 
relationship  to  a  real  structure  in  nature,  and  in  fact  is 
alien  to  the  phenomenon  of  locally  mediated  change.  Rather 
the  "developmental  model"  is  the  aggregate  of  all  its 
submodels  acting  together  fortuitously.  Particular 
developmental  models  must  therefore  be  presented  informally 
as  such  aggregates. 

6.5.1  Modelling  the  Blastoderm  of  Drosophila 

Kauffman  et  al.  (1978)  have  character i zed  the 
blastodermal  stage  of  Drosophila  in  terms  of  a  code  of 
binary  switches  constructed  in  response  to  positive  and 
negative  variations  from  the  uniform  state  of  a  reaction- 
diffusion  system. 


His  model  agrees  with  experimental 
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findings  in  predicting  initial  "segmentation"  of  the 
elongated,  ellipsoidal  egg,  as  sinusoidal  waves  fit  into  its 
length,  followed  as  its  circumference  enters  the  first 
critical  size,  by  division  of  the  egg  into  dorsal  and 
ventral  compartments.  Furthermore  the  codes  thus  assigned 
agree  remarkably  well  with  transdetermi nat ion  data  taken 
from  experiments  such  as  Hadorn' s  (1966),  in  that  the  more 
digits  that  differ,  the  lower  the  tr ansdetermi nat i on  rate. 

Implicit  in  Kauffman's  model  are  certain  assumptions 
concerning  cell  division  in  the  egg.  The  size  of  the  egg 
does  not  increase,  but  nucleation  alters  the  reaction- 
diffusion  parameters,  by  decreasing  diffusion  rates,  in  a 
manner  which  allows  successive  patterns  to  be  supported  in 
the  egg  as  their  wavelengths  decrease.  The  rate  of 
nucleation  is  uniform  throughout  the  egg,  and  this  rate  is 
such  that  odd  patterns  after  the  first  (patterns  asymmetric 
with  respect  to  anterior  end  of  the  egg)  are  missed. 

A  model  has  been  constructed  in  which  these  assumptions 
are  formalized. 

In  the  growth  model  synchronous,  autonomous  cell 
division  is  selected  because  it  (a)  provides  uniform 
expansion,  (b)  results  in  odd  modes  beyond  the  first  being 
missed,  and  (c)  is  mathematically  equivalent  to  variation  of 
system  parameters  as  described  above  (Kauffman  et 
al.  (1978)). 

The  growth  time  step  and  communication  time  step  are 
the  same.  Furthermore,  while  a  general  background  of 


■ 


■ 


197 


unoriented  cell  divisions  is  assumed,  the  orientation  of 
divisions  down  the  egg,  due  to  flows  during  pattern 
formation,  is  modelled  in  order  to  maintain  an  elongated 
shape . 

Since  the  shape  functions  on  ellipses  are  not  yet 
available,  the  blastoderm  is  modelled  here  as  a  rectangular 
structure.  A  code  fragment  "10"  is  added  to  a  cell's 
genetic  code  whenever  a  pattern  arises  in  which  species  one 
is  positive  and  species  two  negative,  and  "01"  when  species 
one  is  negative  and  species  two  positive.  For  purposes  of 
display  these  codes  have  been  translated  into  letters  as 
fol lows : 


A  - 

10 

I  - 

100110 

B  - 

01 

J  - 

100101 

c  - 

1010 

K  - 

011010 

D  - 

1001 

L  - 

011001 

E  - 

0110 

M  - 

010110 

F  - 

0101 

N  - 

010101 

G  - 

101010 

0  - 

010100 

H  - 

101001 

P  - 

010011 

These  codes  translate  directly  to  Kauffman's  setting  of  '  1' 
for  species  one  positive,  '0'  for  species  one  negative. 

The  results  of  simulation  may  be  seen  in  Figure  6.18. 
Note  that  the  width  of  the  system  is  noncritical. 

The  assignment  of  codes  coincides  (under  translation) 
with  Kauffman's  and  therefore  agrees  with  the 
transdetermi nat ion  data.  Segments  are  observable,  though  a 
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FIGURE  6.18a  Segmentation  in  the  blastoderm, 
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FIGURE  6.18b  Assignment  of  code 
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FIGURE  6.18c  Nodal  pattern  two 
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FIGURE  6.18d  Assignment  of  code 
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FIGURE  6 . 1 8e  Nodal  pattern  three 
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mixture  of  two  cell  types  is  seen  near  the  center  of  the 
structure.  This  intermixing  of  cells  is  a  result  of  the 
discrete  representat i on  in  hexagonal  coordinates  and  is 
caused  by  decisions  dependent  on  vertical  lines  being 
imposed  on  the  array.  Any  discrete  representation  must  have 
limitations  in  some  circumstances,  however.  The  results  for 
the  blastoderm,  in  showing  that  segments  do  indeed  arise  as 
determination  occurs  in  response  to  reaction-diffusion  in  an 
elongated  shape,  are  judged  to  be  largely  successful. 

6.5.2  Modelling  the  Chick  Wing 

The  chick  wing  is  a  growing  structure  in  which  cell 
differentiation  into  bone  is  observable.  Initially  a  limb 
bud  appears  and  begins  to  grow  outwards.  During  the 
elongation  phase,  cells  in  the  proximal  central  portion  of 
the  wing  are  observed  "condensing"  into  bone  tissue;  this 
structure  follows  the  process  of  growth  down  the  limb.  At 
the  "elbow"  the  bone  structure  becomes  double,  and  at  the 
wrist  it  separates  into  three  digits. 

Many  models  have  been  proposed  for  the  regulation  of 
these  events.  Wolpert's  (1978)  introduces  the  concept  of 
the  progress  zone  in  which  cells  can  "count"  the  time  they 
spend  in  the  distal  portion  of  the  wing,  and  after  emerging, 
use  this  information  to  differentiate.  Newman  and  Frisch 
(1979)  suggest  a  distal  diffusion  chamber  in  which  the 
pattern  of  cartilage  formation  in  the  plane  perpend i cu 1 ar  to 
the  proximo-distal  axis  of  the  wing  is  specified  by  critical 
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patterns  in  a  reaction-diffusion  system.  Newman  and  Frisch 
suggest  another  mechanism  for  proximo-distal  pattern 
formation.  Wi lby  and  Ede  (1975)  have  developed  a  simulation 
program  in  which  artificially  imposed  periodic  functions  are 
used  as  prepatterns  for  bone  formation.  Moderate  success 
has  been  achieved  in  obtaining  rea 1 i st i c- looki ng  patterns  in 
a  rectangular  coordinate  system.  Wi lby  and  Ede  (1975)  do 
not  here  model  growth  or  i nterce 1 1 u 1 ar  processes  explicitly. 

Newman  and  Frisch's  (1979)  model  with  a  distal 
diffusion  chamber  has  two  flaws.  Firstly,  the  cross-section 
is  modelled  as  rectangular,  thus  being  unable  to  produce  an 
isolated  central  area  of  bone  tissue.  Secondly,  the  model 
has  to  call  upon  a  second  mechanism  to  explain  the  proximo- 
distal  bone  patterns. 

A  model  is  presented  here  which  can  account  for  all  the 
features  of  the  phenomenon  and  can  be  observed  in  action  via 
simulation.  It  is  described  briefly  below: 

(1)  Assume  a  background  rate  of  general  uniform  expansion  of 
the  system. 

(2)  The  initial  configuration  in  the  limb  bud  is  an 
elliptical  cross-section  perpendi cu 1 ar  to  the  proximo-distal 
axis  which  is  of  a  size  critical  for  a  reaction-diffusion 
system  such  as  the  Kauffman  system.  The  shape  function  is 
given  by  the  first  Matthieu  function,  yielding  a  central 
elliptical  core  of  values  negative  in  species  one  and 
positive  in  species  two,  surrounded  by  cells  with  opposite 
signs.  This  configuration  may  be  seen  in  Figure  6.19.  The 
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phase  I 


phase  2 


FIGURE  6.19  Elliptic  cross-sections  in  the  wing 


207 


cross  section  perpendi cu 1 ar  to  the  dor  so- vent ra 1  axis  is 
assumed  to  be  rectangular  and  is  noncritical  in  the  proximo- 
distal  di rect ion . 

(3)  During  general  expansion,  the  elliptical  pattern  is  not 
destroyed,  but  the  proximo-distal  dimension  reaches  critical 
size  for  pattern  1  of  a  mixed  boundary  value  problem.  The 
proximal  boundary  of  the  limb  (continuous  with  body  tissues) 
does  not  obey  a  zero  normal  flux  condition  but  rather 
conforms  to  noncritical  values,  that  is,  to  the  uniform 
stable  state. 

(4)  In  the  proximal  half  of  the  limb  bud,  the  positive 
contribution  from  the  new  pattern  reinforces  the  existing 
pattern  caused  by  the  elliptic  cross-section,  while  in  the 
distal  portion  the  patterns  are  reversed  in  sign. 

(5)  The  rei nforcement  of  the  pattern  means  that  proximal 
tissue  is  held  at  concentration  values  away  from  uniformity 
for  a  prolonged  length  of  time.  Through  either  deprivation 
(causing  repression)  or  saturation  (causing  derepress i on ) 
this  results  in  determination  and  differentiation  of  these 
cells  into  bone  cells,  in  the  central  core,  and  other 
tissue,  in  the  surrounding  cells. 

(6)  The  differentiated  cells  are  no  longer  capable  of 
supporting  the  reaction-diffusion  system,  so  the  system's 
inner  boundary  moves  distal ly. 

(7)  The  system  is  no  longer  in  critical  size  in  the  proximo- 
distal  direction,  so  it  returns  to  uniformity  in  that 
direction  though  still  supporting  the  elliptical  cross- 
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section.  In  returning  to  uniformity  it  causes  cell 
divisions  to  be  oriented  down  the  limb,  thus  causing  it  to 
elongate . 

(8)  The  proximo-distal  dimension  again  becomes  critical  for 
pattern  1,  and  steps  ( 4 ) - ( 8 )  repeat,  until 

(9)  Owing  to  the  background  growth,  the  ellipse  leaves  the 
critical  range  for  the  central  core  pattern  and  enters  the 
critical  range  for  another  pattern  in  which  along  the  major 
axis  of  the  ellipse  ( anter i or -poster ior  axis)  there  are  two 
convex  negative  regions  for  species  one,  while  the  rest  of 
the  ellipse  is  positive  (one  of  the  Matthieu  functions). 
Steps  (3) -(8)  repeat  until 

(10)  During  general  expansion,  the  ellipse  reaches  a  third 
critical  phase  in  which  three  negative  areas  occur  on  the 
major  axis  of  the  ellipse  (another  Matthieu  function). 
Steps  (3) -(8)  repeat  until  growth  is  stopped. 

The  model  has  a  "diffusion  chamber"  or  "progress  zone" 
of  undifferentiated  cells  at  the  distal  end  of  the  limb. 
Surgical  experiments  on  chick  wings  in  which  the  distal 
portion  is  reduced  have  resulted  in  malformed,  randomly 
differentiated  limbs.  The  model  would  predict  such  a  result 
because  elongated  shape  depends  on  pattern  establishment  at 
certain  sizes  and  differentiation  depends  on  pattern  1 
occur ing  at  these  sizes.  So  reducing  the  size  destroys 
these  mechanisms.  The  model  is  sufficient  throughout  the 
development  of  the  wing;  that  is,  all  aspects  of  the  shape 
and  differentiation  phenomenon  can  be  explained  by  this  one 
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mechani sm. 

There  is  an  analogy  to  Wolpert's  clocked  time  in  the 
progress  zone,  since  differentiation  is  provoked  by  a 
certain  "extended  time"  in  a  situation  of  deprivation  or 
saturation.  It  is  felt  however  that  the  present  model  is 
much  closer  to  the  spirit  of  locally  mediated  changes. 
Furthermore  it  explains  geometry  as  well  as  differentiation. 
Limbs  exhibit  elliptic  cross -sect i ons ,  with  more  eccentric 
ellipses  enclosing  more  bone  structures  in  a  proximo-distal 
success i on . 

Prolonged  deprivation  or  saturation  as  a  chemical  basis 
for  differentiation  may  have  significance  in  terms  of  the 
portions  of  DNA  in  which  repeated  sequences  occur.  For 
example,  a  state  of  prolonged  deprivation  of  a  substance, 
which  is  normally  bonded  to  such  an  area  in  multiple  copies, 
might  result  in  evacuation  of  the  substance  to  maintain  the 
standing  pattern.  Removal  of  all  such  bonded  copies,  after 
the  pattern  has  been  in  place  for  a  critical  time  period, 
might  result  in  the  derepression  of  portions  of  DNA  coding 
for  differentiation  into  cartilage.  Such  an  i nterpretat ion 
agrees  with  the  identification  of  repeated  sequences  as 
possible  timing  mechanisms,  which  was  mentioned  in  Chapter 
Three . 

The  simulation  of  the  chick  wing  model  may  be  seen  in 
Figures  6.20-6.22.  Cells  marked  by  "Z"  are  undetermined, 
while  "B"  marks  cartilage  cells  and  "A"  the  surrounding 
tissue.  The  three  phases  of  growth  correspond  to  the  three 
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elliptic  cross  sections 

of 

F i gure  6.19. 

6.6  Conclusion 

The  development 

of 

a  framework  for 

mode  1 

1 i ng  and 

simulating  development 

i  s 

described  in  this 

chapter 

Mode  1 s 

of  cellular  and  i nterce 1 1 u 1 ar  processes  are  simul 

ated,  and 

conditions  under  which 

they  may  be  combined  in 

mode  1 1 i ng 

developmental  systems 

are  discussed. 

These 

processes 

include  cell  division 

> 

the  control  of 

growth 

in  cell 

populations,  intercellular  communication,  pattern  formation, 
determination  and  differentiation.  Within  the  modelling 
structure,  various  options  for  component  models  are 
provided.  The  framework  is  constructed  so  that  at  each 
level  different  components  for  new  phenomena,  or 
at  1 ternat i ves  for  existing  models,  could  be  i ncorporated . 
Application  of  the  modelling  and  simulation  framework  to  two 
developmental  systems,  the  blastoderm  of  Drosophila  and 
chick  wing,  is  described. 


the 
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CHAPTER  SEVEN 


Cone 1 us i ons 


The  primary  goal  of  this  project,  as  stated  in  Chapter 
One,  is  to  provide  a  functioning  and  useful  framework  for 
modelling  in  developmental  genetics.  Along  with 
consideration  of  success  in  achieving  this  goal, 
contributions  to  the  art  of  modelling  and  simulation  and  to 
the  establishment  of  a  mathematical  context  for  development 
can  be  discussed.  Finally,  Ziegler's  theory  of  modelling 
and  simulation  can  be  assessed  in  practical  application  to 
this  type  of  modelling  situation. 

Phenomena  at  various  levels  in  developing  systems  have 
been  modelled  and  simulated  successfully.  The  modelling 
framework  includes  (1)  models  of  cell  division  and  growth 
with  locally  mediated  control  mechanisms,  (2)  models,  both 
numeric  and  analytic,  of  reaction  and  diffusion  as 


mechani sms 


of  pattern  formation,  and  (3)  models  of 


X 
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determination  and  differentiation.  Two  developmental 
systems,  employing  a  combination  of  these  modelling  options, 
are  presented.  These  are  segmentation  in  the  blastoderm  of 
Drosophila  and  bone  condensation  and  growth  in  the  chick 
wing . 


Some 


observations  should  be  made  about 


the 


contributions  of  these  developmental  models.  Firstly,  the 
various  phenomena  and  processes  of  development  have  been 
modelled  and  simulated  in  interaction,  as  is  possible  only 
through  a  framework  in  which  different  models  can  be 
combined.  Secondly,  all  models  in  the  framework  adhere 
strictly  to  the  requirement  that  changes  be  locally 
mediated.  Thus,  for  example,  departures  from  initially 
uniform  concentrat ion  configurations  are  not  induced  by  an 
imposed  polarization  of  some  kind,  but  by  a  local  process 
(reaction-diffusion)  considered  in  an  i nterce 1 1 u 1 ar  context 
(the  morphogenetic  field).  Similarly,  cells  experience 
genetic  changes  individually  in  response  to  critical  values 
of  local  influences.  Fundamental  to  the  ability  of  the 
framework  to  model  interacting  processes,  and  locally 
mediated  change,  is  the  provision  of  modelling  structures 
which  allow  sufficient  information  to  be  maintained  for  each 
cell,  and  permit  modelling  to  be  done  at  higher  levels  of 
the  developmental  system. 

At  each  level  of  modelling  there  are  opportunities  for 
further  options  to  be  developed.  Different  models  of 


growth,  both  intracel lular  and  i nterce 1 1 u 1 ar , 


could 


be 
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formulated.  Intercellular  communication  other  than  by 
reaction  and  diffusion  could  be  studied  (if  other  stable 
mechanisms  were  found).  Determination  and  differentiation 
in  a  variety  of  one  and  two  dimensional  systems  could  be 
modelled;  for  example,  experiments  in  the  ligation  of  insect 


eggs  could  be 

per  formed 

( Kauffman  et  a  1 . , 

1978) 

It  i  s  a 

measure  of  the 

power  and 

flexibility  of 

the 

mode  1 1 i ng 

framework  that  such  extensions  are  immediately  feasible. 

The  modelling  framework  is  significantly  limited,  in 
its  present  state  of  development,  by  difficulties  in 
studying  the  establishment  of  patterns  in  non- rectangu 1 ar 
domains.  On  the  one  hand,  numerically  integrated  solutions 
in  such  domains  are  expensive  (though  feasible)  to  obtain. 
On  the  other  hand,  the  shape  functions  with  which  analytic 
solutions  can  be  calulated  may  be  unknown  (as  in  irregular 
domains)  or  difficult  to  use  (as  in  circular  or  elliptic 
domai ns ) . 

It  may  be  noted,  however,  that  if  numerical  integration 
methods  are  used,  they  are  expensive  only  during  critical 
phases,  when  the  system  is  either  leaving  or  returning  to 
the  uniform,  steady  state  configuration.  During  those 
periods  of  time  when  the  system  is  noncritical,  convergence 
to  the  steady  state  solution  is  immediate  at  each  time  step. 
Furthermore  patterns  such  as  those  made  up  of  Bessel  and 
Mathieu  functions  may  be  obtainable  analytically  by  using 
numerical  techniques  to  generate  polynomial  approximations 
in  their  respective  domains.  The  scope  of  the  modelling 
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framework  would  be  increased  considerably  by  the  addition  of 
such  techniques.  Modelling  embryonic  development,  imaginal 
discs  in  Drosophila,  and  other  systems  decomposible  into 
(approximately)  regular  domains,  would  then  become  feasible 
in  the  framework. 

The  data  structures  of  the  modelling  framework  could  be 
employed  fruitfully  in  other  biological  simulations  where 
the  capacity  for  dynamic,  locally  mediated  changes  in  cell 
space  structures  is  required.  Such  areas  include  studies  of 
cell  movements,  cell  sorting  models,  and  brain  models.  It 
may  be  noted  with  regard  to  the  last  that  reaction-diffusion 
systems  are  being  investigated  in  modelling  memory 
(Tuckwell,  1979).  Experiments  with  cancer  and  clonal 
studies  of  mutated  cells  can  also  be  studied  in  this 
f r amework . 

The  framework  is  thus  able  to  accommodate  a  resonable 
variety  of  phenomena  in  the  biological  domain.  This 
capacity  is  felt  to  be  an  indication  of  structural  validity 
in  the  many  component  models  of  the  framework  and  in  the 
ways  in  which  they  can  be  combined. 

Simulation  of  the  models  in  the  framework  makes  special 
contributions  to  developmental  modelling.  Firstly,  in  order 
for  behavioral  data  to  be  generated  by  simulation,  every 
aspect  of  that  simulation,  and  of  the  models  it  represents, 
must  actually  perform  as  expected  operationally.  Thus  in 
the  model  of  Section  6.5.2  pre-pat terni ng  and 
differentiation  in  the  chick  wing  can  be  observed  in  action. 
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A  faulty  model  such  as  that  of  Newman  and  Frisch  (1979) 
could  not  generate  such  behavior.  In  contrast,  the  model  of 
Kauffman  et  al.  (1978)  of  the  blastoderm  does  generate  its 
predicted  behavior.  Simulation,  when  performed  with  regard 
to  validity  of  the  iterative  representat ion  of  a  model, 
fulfils  the  function  of  ensuring  that  the  model  is  a  valid 
generator  of  the  desired  behavior.  Secondly,  through 
simulation  an  inferred  process  Known  as  "reaction  and 
diffusion"  has  been  demonstrated  to  be  a  viable  tool  in 
pattern  formation.  Thirdly,  the  vital  importance  of  the 
relative  speed  with  which  growth  takes  place  and  reaction- 
diffusion  operates  has  been  illustrated  by  means  of 
simulation.  Which  patterns  appear  in  a  system,  and 
therefore  what  structures  of  differentiated  and  determined 
cells  will  occur,  depend  on  this  timing.  It  is  suggested 
that  the  interaction  of  these  phenomena  (growth  and 
communication),  in  combination  with  viability  tests  in  the 
process  of  natural  selection,  may  play  a  part  in  the 
evolution  of  species. 

In  view  of  the  above  contributions,  it  would  appear 
that  a  useful,  functioning  framework  for  developmental 
modelling  has  been  largely  achieved.  Contributions  to  the 
art  of  modelling  and  simulation  may  now  be  discussed. 

The  modelling  framework  is  mu  1 1 i 1  eve  1 1 ed :  modelling  is 
possible  at  the  molecular  level  (reaction  and  diffusion  of 
chemical  species,  for  example),  at  the  cellular  level  (cell 
division,  for  example),  and  at  the  level  of  the 
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developmental  system  (the  chick  wing,  for  example).  The 
framework  is  also  of  mixed  type:  state  variables  may  be 
discrete  (clocks),  binary  digit  strings  (genetic  code),  or 
real  numbers,  (values  in  a  continuous  field  such  as  local 
concentrations).  The  framework  forms  a  "model  pool"  of  the 
type  envisaged  by  Ziegler  (1976),  and  it  has  been 
demonstrated  how  selected  members  of  the  pool  act  together 
in  obedience  to  conditions  of  validity  established  with 
respect  to  their  common  time  base. 

The  modelling  framework,  therefore,  is  a  functional 
example  of  aspects  of  modelling  and  simulation  of  current 
interest.  The  techniques  developed  to  deal  with  dynamic 
change  in  cell  structures  and  locally  mediated  processes 
(local  transition  functions)  may  prove  valuable  in  any  area 
of  computing  science  where  the  cell  space  formalism  is  used. 
In  summary,  it  is  hoped  that  contributions  have  been  made  to 
the  art  of  modelling  and  simulation. 

The  mathematical  context  which  has  been  presented  for 
development  is  a  systems  theoretic  formalism  operating  in  a 
cell  space.  Through  its  components  (cells  and  lists)  it 
proves  to  be  sufficient  in  representing  the  structures  of 
the  chosen  experimental  frame.  Through  its  transition 
functions  the  formalism  appears  capable  of  representing 
dynamic  processes  in  these  structures.  Such  processes 
employ  those  mathematical  tools,  pertaining  to  various 
levels  of  the  system,  which  are  identified  in  Chapter  One  as 
being  necessary  components  of  developmental  modelling.  For 
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example,  differential  equations  are  employed  at  the 
molecular  level,  statistical  techniques  are  used  in 
modelling  cell  growth  in  a  population,  and  equilibrium 
phenomena  are  employed  in  morphogenetic  fields. 

The  formalism  proves  to  be  capable  of  representing 
dynamic  changes  in  structures,  that  is  in  representing  the 
interaction  of  the  dynamic  and  static  elements  of  the 
modelling  domain.  An  example  of  change  in  static  elements 
or  structure  is  given  by  growth,  which  may  be  controlled  by 
dynamically  changing,  locally  mediated  processes.  An 

example  of  change  in  dynamic  elements  is  given  by 
determination,  in  which  the  genetic  code  changes  in  response 
to  critical  values  of  some  local  phenomena,  resulting  in 
changed  act i vi ty . 

Simulation  of  the  system  specifications  forms  a  vital 
part  of  the  mathematical  context  which  has  been  employed. 
Dynamic  aspects  of  models  without  simulations  can  only  be 
appreciated  by  means  of  inference  from  the  transition 
function  specifications,  whereas  simulation  actually 
generates  behavior  for  the  model.  The  transition  functions 
can  thus  be  observed  in  action,  and  since  they  are  meant  to 
describe  action  this  is  the  appropriate  technique. 

Questions  of  validity  must  of  course  be  addressed;  this  may 
be  done  directly  since  Ziegler's  (1976)  theory  of  modelling 
and  simulation  operates  in  the  same  formalism.  In  many 
respects,  therefore,  the  mathematical  context  chosen  for 
modelling  the  developmental  system  proves  to  be 
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sat i sf actory . 

Ziegler's  theory  of  modelling  and  simulation  (Ziegler, 
1976)  has  been  used  throughout  in  the  development  of  the 
modelling  framework  and  in  its  exposition.  The  formal 
specifications  of  models  in  Chapter  Six  proved  to  be 
cumbersome  at  times,  yet  this  is  justifiable  on  two  counts. 
Firstly,  sets  of  component  models  of  the  "modelling  pool" 
can  be  selected  with  regard  to  formally  established  criteria 
for  vali'dity  (as,  for  example,  in  the  selection  of  a  growth 
and  signalling  model  pair).  Secondly,  changes  in  structure 
and  local  functioning  are  describable  in  systems  theoretic 
terms  in  a  manner  which  relates  them  to  global  structure; 
this  perspective  could  not  be  gained  if  they  were  to  be 
described  simply  as  individual  algorithms. 

Morphisms  between  models  at  the  iterative  level  have 
not  proved  useful.  It  is  felt  that  they  might  become 

applicable  if  the  modelling  framework  were  to  be  extended. 
The  primary  thrust  of  a  project  such  as  this  is  in 
establishing  the  different  levels  of  the  modelling  framework 
(a  "vertical"  phase  of  development).  Later  expansion,  by 
means  of  other  models,  and  lumped  models  of  existing  models, 
would  occur  at  each  level  (a  "horizontal"  phase  of 
development).  Relationships  between  models  on  a  particular 
level  could  then  be  investigated,  by  means  of  canonical 
decompositions  of  input  segments  and  the  relevant  morphisms. 

Ziegler's  theory  has  proved  invaluable  as  a  means  of 
keeping  perspective  in  the  construction  of  a  modelling 
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framework  of  this  complexity  and  scope.  Organized  progress 
from  the  real  system  to  simulation  experiments,  as  followed 
in  the  thesis  and  in  the  development  of  the  modelling 
framework,  with  attention  paid  to  questions  of  validity  at 
each  stage,  has  resulted  in  a  functioning,  effective 
modelling  framework. 

In  summary,  it  may  be  stated  that  many  of  the  goals  of 
the  project  have  been  accomplished.  A  useful  framework  for 
modelling  and  simulation  in  developmental  genetics  has  been 
accomplished,  with,  it  may  be  hoped,  benefit  to  both  biology 
and  computing  science. 
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